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Abstract 


While the performance of flight simulator motion system hardware has advanced 
substantially, the development of the motion cueing algorithm, the software that 
transforms simulated aircraft dynamics into realizable motion commands, has not kept 
pace. Prior research identified viable features from two algorithms: the nonlinear 
“adaptive algorithm”, and the “optimal algorithm” that incorporates human vestibular 
models. A novel approach to motion cueing, the “nonlinear algorithm” is introduced that 
combines features from both approaches. This algorithm is formulated by optimal 
control, and incorporates a new integrated perception model that includes both visual and 
vestibular sensation and the interaction between the stimuli. Using a time-varying 
control law, the matrix Riccati equation is updated in real time by a neurocomputing 
approach. 

Preliminary pilot testing resulted in the optimal algorithm incorporating a new 
otolith model, producing improved motion cues. The nonlinear algorithm vertical mode 
produced a motion cue with a time-varying washout, sustaining small cues for longer 
durations and washing out large cues more quickly compared to the optimal algorithm. 
The inclusion of the integrated perception model improved the responses to longitudinal 
and lateral cues. False cues observed with the NASA adaptive algorithm were absent. 
The neurocomputing approach was crucial in that the number of presentations of an input 
vector could be reduced to meet the real time requirement without degrading the quality 


iii 


of the motion cues. 



The new cueing algorithms are implemented on the NASA Langley Visual 
Motion Simulator (VMS), and will ultimately be implemented on the new Cockpit 
Motion Facility (CMF) currently being erected at NASA Langley. 
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1. Introduction 


1.1. Effect of Simulator Motion on Human Performance 

The objective of a motion system, when used in conjunction with a visual system, 
is to stimulate the pilot so that he or she can perceive the required motion and force 
information (i.e., cues) necessary to fly the simulator within the same performance and 
control activity as the actual aircraft. An example of a motion system is the six-degree- 
of-freedom hexapod shown in Figure 1.1. 



Figure 1.1. Six-Degree-of-Freedom Hexapod Motion System. Delft University, The 
Netherlands. 

Buckingham [1] reported that the inclusion of motion cues allows the pilot to 
become aware of the aircraft response before visual cues are detected, noting that without 
motion cues, the pilot’s perception of motion is degraded and the aircraft feels slower in 
responding. Buckingham noted that in extreme cases, the pilot might be unable to control 
the aircraft when the absence of motion cues introduces a 90-degree phase lag into the 
control loop. Buckingham cited one case in which the motion system was disabled while 
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unbeknown to the pilot; the pilot described the simulation as inferior to the previous 
configuration with motion, as it took longer to respond to control inputs. 

Gundry [2] reported that Douvillier, et al., Fedderson, Mathney, Perry and Naish, 
and Tremblay, et al. observed that when motion cues were provided, there was an 
increase in high-frequency, low-amplitude control movements that appeared more like 
movements observed during flight as compared to activity in a fixed-base simulator. 
Gundry noted that Perry and Naish compared pilot control activity of both fixed-base and 
moving-base simulation of a flight through heavy turbulence, and observed a 
considerable reduction in the simulated aircraft roll angle with motion present. The 
presence of motion produced pilot responses with more rapid and accurate control. These 
results show that when an aircraft is subjected to turbulence in flight, the pilot uses roll 
and pitch motion as information to correct the aircraft attitude. It was observed that 
platform motion in response to external disturbances and maneuvering allowed the 
operator to control the simulator using sensory cues similar to those used in flight. 

Gundry also reported that Dinsdale, Meiry, Shirley, and Stapleford investigated 
the effects of roll motion upon compensatory tracking error. In these investigations, the 
presence of motion was observed to reduce the phase lag of the simulated aircraft roll 
angle relative to the command input, increase mid-frequency gain and crossover 
frequency, and reduce the size of the remnant. This showed that the presence of roll 
motion cues provide the operator with lead information that is used to track the 
disturbance input more accurately, especially at frequencies greater than 0.5 Hz (as noted 
by Shirley). 


2 



Scanlon [3] conducted a piloted simulation study on the NASA Langley Visual 
Motion System to determine the effects of motion cues during the performance of 
complex curved approach and landing tasks in the signal environment of the Microwave 
Landing System (MLS). Comparisons of pilot tracking performance and workload were 
made on approach tasks of low, medium, and high complexity conducted with and 
without motion, with and without turbulence, and with three different wind models. With 
motion cues, smaller lateral tracking errors resulted for the most complex approach in the 
presence of wind and turbulence. The effect of motion was insignificant for lateral 
tracking errors for low and medium complexity approaches, and for vertical tracking 
error for all levels of complexity. Motion cues, most noticeably with turbulence, yielded 
a higher physical workload as measured by pilot control activity, with higher column and 
wheel input rates measured for all levels of task complexity. All pilots indicated a 
preference for motion over no motion, commenting that flying was easier and more 
realistic with the addition of motion. 

Schroeder [4] conducted an evaluation on the NASA Ames Vertical Motion 
System with experienced test pilots performing single-axis vertical and directional (yaw) 
maneuvers of a hovering helicopter with varying degrees of fidelity in the motion cueing 
algorithm, i.e., from nearly full motion to fixed-base. For the vertical maneuver full- 
motion case, Schroeder reported that “well-damped, accurate bob-ups are achieved with 
the vertical velocity staying within 10 ft/s”, but for the fixed-base case, the pilot had to 
adjust his compensation with the remaining cues, recovering over time, but taking longer 
to achieve final repositioning. Schroeder noted that the pilots were “stunned” by the total 
loss of motion, reporting that motion cues were “certainly perceived” by all pilots for all 
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test conditions with motion cues, but not with the fixed-base configuration. Schroeder 
commented, “Until the value of motion was demonstrated, pilot subjective impression 
was that the vertical task was primarily visual”. Schroeder reported for the directional 
maneuver, no performance degradation was noticed for the fixed-base case, noting that 
visual yaw cues, depending on the visual scene, may be very compelling in inducing 
motion perception, possibly overwhelming the yaw motion stimulus. 

Hall [5] noted that platform motion remains the only currently available 
technology that can provide motion cueing of both direction and magnitude without 
requiring additional learning, because the pilot’s proprioceptive sensors are stimulated in 
the short term in the same manner as in flight. The presence of motion will allow the 
pilot to achieve a task closer to that seen in the aircraft since he uses a similar set of 
sensory cues, especially when forced to operate in a high gain manner. Hall then 
summarizes that motion becomes less important when the vehicle is easy to fly, the task 
can be performed with low pilot workload and gain, and disturbance motion is either 
absent or does not require corrective action. Platform motion becomes increasingly 
important as task difficulty and pilot control gain increase, and are essential in the 
absence of good, wide field of view visual cues (e.g., flying in clouds, at night), and 
necessary for high gain tasks even with strong visual cues. Hall concluded that motion 
cueing is essential when a pilot must either react quickly in response to an unexpected 
disturbance, or when the pilot must control a vehicle with low stability. 

1.2. Vehicle Simulation Structure 

The vehicle simulation structure for a motion system is shown in Figure 1.2. The 
operator control inputs drive a mathematical model of the vehicle dynamics, generating 
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the vehicle states. Passing the vehicle states through the motion cueing algorithm 
produces the desired motion cues and platform states. The desired platform states are 
then transformed from degree-of-freedom space to actuator space, generating the realized 
commands to the six actuators. The actuator motion commands serve as input to the 
platform dynamics, resulting in the actual simulator motion. 



Figure 1.2. Vehicle Simulation Structure. 

The motion cueing algorithm generates the desired motion cues that are 
constrained within the physical limits of the motion system. Figure 1.3 shows a typical 
motion cueing algorithm implementation. Vehicle states are transformed from a body 
reference frame to an inertial reference frame. Scaling and limiting the vehicle states 
reduces the magnitude of the motion cues. The duration of the cues are limited by the 
physical dimensions of the motion system. A method to overcome this limitation is a 
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technique known as “washout”. Washout involves returning the platform state to a 
neutral position following the initial, or “onset” portion of a motion cue, thus “washing 
out” the resulting cue at levels below the pilot’s perceptual threshold. This is 
accomplished by passing the vehicle state through a high-pass filter, removing long- 
duration (low-frequency) motion components. Figure 1.4 shows the response of a high- 
pass washout filter to an acceleration ramp to step input. 


Pa 



Figure 1.3. Motion Cueing Algorithm Implementation. 

The otolith organs in the human vestibular system sense both acceleration and 
tilting of the pilot’s head with respect to the gravity vector. Since the otoliths cannot 
discriminate between acceleration and tilt, this phenomenon, known as tilt coordination, 
can be used to advantage in motion simulation. For long-term specific force simulation, 
tilting the motion platform at a rate below the pilot’s perceptual threshold augments the 
short-duration acceleration cues produced by high-pass washout filters. This additional 
cue results from passing the vehicle acceleration through a low-pass filter to produce the 
desired long-duration tilt cue. Tilt coordination is implemented in a motion cueing 
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algorithm by adding additional cross-feed channels with low-pass filters in the 
longitudinal (pitch/surge) and lateral (roll/sway) modes that produce the additional 
rotational cues as shown in Figure 1.3. For this reason four separate modes are 
implemented in a motion cueing algorithm: longitudinal, lateral, yaw, and heave. 
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Figure 1.4. Response of a High-Pass Washout Filter to a Ramp to Step Input. 


1.3. Scope of Research 

In recent years, the performance of the hardware used to create the sensation of 
motion in flight simulators has improved substantially. However, development of the 
motion cueing algorithm, the software that transforms the simulated aircraft dynamics 
into realizable commands to the motion system hardware, has not kept pace with the 
hardware development. Wu and Cardullo [6] reported that early approaches to motion 
cueing using simple (first- and second-order) linear washout filters, for which the ratio of 
onset to washout duration was fixed, resulted in poor motion cues. This was a 
consequence of the ratios of onset to washout duration and magnitude being fixed, 


7 





thereby limiting the duration of low-magnitude cues to that of the maximum cue. In 
addition, most existing algorithms are oriented towards minimizing the state error 
between the aircraft and simulator rather than the perceptual error between the aircraft 
and simulator pilot. Wu and Cardullo [6] identified the two most viable approaches to 
motion cueing, the nonlinear “adaptive algorithm”, for which the ratios of onset to 
washout duration and magnitude vary with time, and the linear, human-centered “optimal 
algorithm”. 

The coordinated adaptive washout algorithm, or “adaptive algorithm” was 
developed at NASA [7]. The objective of this algorithm is to adjust the motion platform 
response based upon its current motion states by adjusting filter gains through a process 
of minimizing a cost function in real time. The cost function is minimized by 
continuously adjusting a set of adaptive parameters by the method of steepest descent. 
This technique has at its basis the minimization of state error between the aircraft and 
simulator. This algorithm is described in further detail in Section 2.7. 

The “optimal algorithm” was developed by Sivan, et al. [8], and later 
implemented at the University of Toronto Institute of Aerospace Studies (UTIAS) [9, 10]. 
This algorithm uses higher-order linear filters that are developed, prior to real time 
implementation, using optimal control methods. This method incorporates a 
mathematical model of the human vestibular system, constraining the pilot sensation 
error between the simulated aircraft and motion platform dynamics. Wu and Cardullo [6] 
reported that the optimal algorithm showed the most potential for future research, 
although the time- varying feature of the adaptive algorithm was also desirable. 
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A primary component in motion simulator design is the determination of the 
motion information that is relevant to the task and has an impact on human performance. 
This requires knowledge of human motion perception that, when integrated in the cueing 
algorithm development can provide the most necessary and beneficial motion cues. To 
that end, an integral part of this research involved the modeling of the human vestibular 
and perceptual systems. Literature studies in motion sensation and the vestibular system 
have been conducted to develop vestibular system sensation models that are most 
consistent with both experimental and theoretical analyses. New models of both the 
semicircular canals and the otoliths are proposed. Literature studies of the characteristics 
of visually induced motion sensation and the visual-vestibular interaction have also been 
conducted. A new integrated human perception model is proposed that includes both 
visual and vestibular sensation and incorporates the interaction between the stimuli. 

The new vestibular models are incorporated in an improved development of the 
linear optimal algorithm. The development of this algorithm is presented along with 
results that demonstrate the effects of implementing the new vestibular models. The 
nonlinear algorithm is a novel approach to motion cueing that combines features of the 
nonlinear adaptive and linear optimal algorithms. This algorithm incorporates the human 
vestibular models along with the new integrated human perception model. The algorithm 
is formulated as an optimal control problem with a nonlinear control law, resulting in a 
set of nonlinear cueing filters that are adjusted in real time based on the motion platform 
states. A neurocomputing approach to solve the matrix Riccati equation in real time is 
discussed. Responses to single degree-of-freedom aircraft inputs for the nonlinear 
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algorithm are presented in comparison with the NASA adaptive algorithm and the 
optimal algorithm. 
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2. Background Information 

2.1. NASA Langley Visual Motion Simulator (VMS) 

The NASA Langley Visual Motion Simulator (VMS), shown in Figure 2.1, is a 
general-purpose flight simulator consisting of a two-crewmember cockpit mounted on a 
60-inch stroke six-degree-of- freedom synergistic motion base [11], [12]. 



Figure 2.1. NASA Langley Visual Motion Simulator (VMS). NASA Langley 
Research Center, Hampton, Virginia. 

Motion cues are provided in the simulator by the extension or retraction of the six 
hydraulic actuators of the motion base relative to the simulator neutral position. The 
NASA adaptive algorithm and the new optimal and nonlinear algorithms were used to 
drive the motion base during the tuning of the new algorithms and the piloted test 
evaluation. 

The cockpit of the VMS, shown in Figure 2.2, is designed to accommodate a 
generic transport aircraft configuration on the left side and a generic fighter or rotorcraft 
configuration on the right side. Both sides of the cockpit are outfitted with three heads- 
down CRT displays (primary flight display, navigation/map display, and engine display), 
a number of small standard electromechanical circular instruments and a landing gear 
handle mounted in the instrument panel. The left side contains a two-axis side stick 
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control loader, and the right side contains a control loaded two-axis center stick. Both 
sides contain control loaded rudder systems. The center aisle stand is outfitted with a 
control display unit, a four-lever throttle quadrant, a flap handle, a speed brake handle, 
and a slats handle. The cockpit is outfitted with four collimated window display systems 
to provide an out-the-window visual scene. During the piloted evaluations, the test 
subject flew from the left seat, while an observer/test conductor rode in the right seat. 



Figure 2.2. Visual Motion Simulator Cockpit. NASA Langley Research Center, 
Hampton, Virginia. 

The simulator includes a high fidelity, highly nonlinear mathematical model of a 
Boeing 757-200 aircraft, complete with landing gear dynamics, gust and wind models, 
flight management systems, and flight control computer systems. For this study, the test 
subjects flew the simulated aircraft in the manual control mode (without the autopilot), 
and with manual throttle control (without the auto throttle). 

The out-the window visual scene is driven by an Evans and Sutherland ESIG 
3000/GT computer generated image system. The visual database represented the 
Dallas/Fort Worth airport and its surrounding terrain. The study utilized runways 18L 
and 18R for approach maneuvers and runway 18R for takeoff maneuvers. The runways 
were equipped with approach lights, precision approach path indicator lights, runway 
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markings, and signage. The database included all runways and taxiways, and all airport 
structures and buildings. All tests were conducted in a daylight environment with full 
visibility. 

2.2. Reference Frames 

A series of reference frames are used in the definition of the motion cueing 
algorithms. These reference frames are defined below and are shown in Figure 2.3. 

2.2.1. Aircraft Center of Gravity 

The aircraft center of gravity reference frame Fr CG has its origin at the center of 
gravity of the aircraft. Frame Fr CG has an orientation for X CG , Y CG , and Z CG that is parallel 
to reference frames Fr s and Fr A . 

2.2.2. Simulator 

The simulator reference frame Fr s has its origin at the centroid of the simulator 
payload platform, i.e. the centroid of the upper bearing attachment points. The origin is 
fixed with respect to the simulator payload platform. X s points forward and Z s points 
downward with respect to the simulator cockpit, and Y s points toward the pilot's right 
hand side. The x-y plane is parallel to the floor of the cockpit. 

2.2.3. Aircraft 

The aircraft reference frame Fr A has its origin at the same relative cockpit location 
as the simulator reference frame Fr s . Fr A has the same orientation for X A , Y A , and Z A 
with respect to the cockpit as the simulator frame Fr s . 
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2.2.4. Inertial 


The inertial reference frame Fri is earth-fixed with Z r aligned with the gravity 
vector g. Its origin is located at the center of the fixed platform motion base. Xi points 
forward and Yi points to the right hand side with respect to the simulator pilot. 

2.2.5. Reference Frame Locations 

In Figure 2.3 are four vectors that define the relative location of the reference 
frames. Ri defines the location of Fr s with respect to Fri. R s defines the location of Frps 
with respect to Fr s . Similarly, R A defines the location of Fr PA with respect to Fr, v where 
R a = R s . R cg defines the location of Fr A with respect to Fr CG . 



Figure 2.3. Reference Frame Locations. Adapted from Wu [13]. 
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2.3. Coordinate Transformations 

The orientation between the body-fixed simulator reference frame Fr s and the 
inertial reference frame Fr, can be specified by three Euler angles: p = [(/) 9 !//]' that 
define a sequence of rotations that carry Fr s into Fiy A vector V expressed in the two 
frames can be related by the transformation matrix L IS (Fp to Fr s ) or L SI (Fr s to Fr,), with 
V s = L IS V‘ and V 1 = L SI V S , where L IS = , and 

cos9cosy/ sin^sin^cos^-cos^sin^ cos^sin^cos^ + sin^sin^r 
L SI = cos#sin^/ r sin^sin^sin^ + cos^cos^ cos0sin#sin^/'-sin0cos^ . (2.1) 

-sin# sin^cos# cos^cos# 

The angular velocity of Fr s with respect to Fr, can be related to the Euler angle 

rates p by the following expression. Let cu* represent the components of this angular 

velocity in frame Fr s , then P = T s to* , where 

1 sin 0 tan 0 cos 0 tan 9 

T s = 0 cos^ -sin^ , (2.2) 

0 sin 0 sect? cos 0 sec# 

and to* = T s _1 p , where 

1 0 -sin# 

T s 1 = 0 cos^ sin^cos# . (2.3) 

0 -sin^ cos^cos# 

Note that in this example, the body-fixed aircraft reference frame Fr A can replace the 
body-fixed simulator reference frame Fr s . 
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2.4. Actuator Geometry 

The geometry of a six-degree-of-freedom synergistic motion system is given in 
Figure 2.4. The relevant vectors relating the locations of the upper and lower bearings of 
the j-th actuator are shown below in Figure 2.5. 



Figure 2.4. Geometry of a Six-Degree-of-Freedom Motion System. Adapted from 
Wu [13]. 
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Figure 2.5. Vectors for the j-th Actuator. 


In Figure 2.5, O s and Oi are the centroids of the motion platform and fixed 
platform respectively, and are also respectively the origins for Fr s and Fiy It can be seen 
that the relation among those vectors is 

R I+ A;= Rj =B; + Z r (2.4) 

The actuator length vector can then be found from 

h= A )+ R i~K ( 2 - 5 ) 


The expression of /. in the inertial reference frame Fri is desired: 


= L si Aj + Rj-B], 


(2.6) 


where A? are the coordinates of the upper bearing attachment point of the j-th actuator in 


Frs and B‘ are the coordinates of the lower bearing attachment point of the j-th actuator 
in Ff|. The actuator extensions can then be found from 


a/‘=/;m-/;(o) 

= (M‘) + ( R . (*) - R . (°)i 


(2.7) 
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Usually the actuator extension is computed from a neutral platform position, where 


L S i(0) = I ( I is the identity matrix), and R , = 0 . therefore A/. 1 = L SI Aj + AR, . 

2.5. Nonlinear Input Scaling 

Limiting and scaling are applied to both aircraft translational input signals a^ and 
rotational input signals . Limiting and scaling modify the amplitude of the input 
uniformly across all frequencies. Limiting is a nonlinear process that clips the signal so 
that it is limited to be less than a given magnitude. Limiting and scaling can be used to 
reduce the motion response of a flight simulator. A third-order polynomial scaling was 
developed [13] and has been implemented in the new simulator motion cueing 
algorithms. 

When the magnitude of input to the simulator motion system is small, the gain is 
desired to be relatively high, or the output may be below the pilot’s perception threshold. 
When the magnitude of input is high, the gain is desired to be relatively low or the 
simulator may attempt to go beyond its hardware limits. Let us define the input as x and 
the output as y. Now define x max as the expected maximum input and y max as the 
maximum output, and ,v 0 and s 1 as the slopes at x = 0 and x = x max respectively. Four 
desired characteristics for the nonlinear scaling are expressed as: 


(1) x = 0 => y = 0, 

(2) ^ ~ "^max y ~ y max ’ 



A third-order polynomial is then employed to provide functions with all the 
desired characteristics. This polynomial will be of the form 
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where 


y = c 3 x 3 + c 2 x 1 + c,x + c 0 , 


(2.8) 


c 0 = 0, 

Cl = V 

C 2 ~ Tnax (^Jmax _ ^ S Q X max ~ ‘h^max )’ 

C 3 = '^max ( Vmax “ 2 )Vax + ^max )• 

One example of this polynomial gain is shown in Figure 2.6, with parameters set as 

*max= 10 > ymax = 6 > 5 0 = 1 A 5, = 0.1. 
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Figure 2.6. Nonlinear Input Scaling. 

2.6. Specific Force at the Pilot’s Head 

The purpose of the motion cueing algorithm is to create a specific force vector 
and an angular velocity vector at the pilot's location in the simulator that approximates 
the stimulus that the pilot would experience in an actual aircraft. The relation between 
the specific force acting on the simulator pilot and the specific force at the origin of the 
simulator reference frame can be found from 


Nonlinear Scaling 
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(2.9) 


fps a ps 8 

=a|? + tb^ x R s + co^ x (cd; 5 x R s ) - § S 
=f s S + CD* x R s + ®s x (c>s x R s ). 

Both fp S and co|! are used to compute sensed responses using the vestibular models 

discussed in Chapter 3. Similar expressions can be obtained for the specific force and 
angular velocity at the aircraft pilot’s head. 

2.7. Coordinated Adaptive Washout Algorithm 

The intent of the NASA adaptive algorithm [7] is to adjust the response of the 
simulator washout filters in real time according to the current state of the simulator. The 
block diagram for this algorithm is shown in Figure 2.7. There are separate filtering 
channels for the translational and rotational degrees of freedom with a cross-feed path to 
provide the steady-state tilt coordination cues. 


Ps 



Figure 2.7. Coordinated Adaptive Washout (NASA Adaptive) Algorithm. 
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The aircraft acceleration vector a A is first transformed from the center of gravity 
of the aircraft to the motion base centroid. After nonlinear scaling and limiting, the 
gravity vector is subtracted to produce a simulator frame specific force vector. The 
simulator specific force is transformed from the simulator frame Fr s into the inertial 
frame Fr h resulting in the inertial specific force command f' . The specific force 
command f{ is passed through a translational channel with a time- varying gain A to 
produce a simulator translational acceleration command S 1 . This acceleration is 
integrated to produce the velocity S 1 , which is then integrated to produce the simulator 
translational position command S 1 . Both the velocity and position commands are 
employed as feedback. 

The aircraft angular velocity vector co A is limited and scaled similar to the 
translational channel, with the resulting vector being transformed to the Euler angular 
rate vector P A . This vector is passed through the rotational channel with a time-varying 
gain 8 to produce the vector P SR . The tilt coordination rate P ST is formed from the 
acceleration a A being passed through the cross-feed channel with a fixed gain y. The 
summation of P ST and P SR will yield P s , which is then integrated to generate P s , the 
simulator angular position command. 

Lsi and T§ are formed by Eqs. (2.1) and (2.2). The simulator translational 
position S 1 and the angular position P s are used to transform the simulator motion from 

degree-of-freedom space to actuator space as given in Eqs. (2.6) and (2.7), generating the 
actuator commands required to achieve the desired platform motion. 
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The control law for the longitudinal mode is given by the following expressions: 


K = KfL-d x S[-eX 

4 = yJL + 44 . 


( 2 . 10 ) 


where d x , e x , and y x are fixed parameters, and A x and S x are the time-varying 

parameters that are continuously adjusted in an attempt to minimize the instantaneous 
value of the cost function. The cost function is defined as 


J, =\(fL -4'f + ^(4 -4f + b Xs',) 2 +^(s!) 2 . 


( 2 . 11 ) 


where W x , b x , and C x are constant weights that penalize the difference in response 

between the aircraft and simulator, as well as restraining the translational velocity and 
displacement in the simulator. 

The time-varying parameters A x and 8 X are adjusted by steepest descent as given 


by 




4 ) 


( 2 . 12 ) 


where K x , K a , K s , and K jS are constants. The first right-hand side term of each 

equation defines the change of the time-varying parameter is to be toward a minimum, 
and together with the second term defines the rate of change. The second term also 
restrains the deviation of either A x or S x respectively from their original values. 
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3. Vestibular System Modeling 

This chapter discusses the development of vestibular system sensation models that 
are most consistent with both experimental and theoretical analyses and can be readily 
implemented into a motion cueing algorithm. These results are based on the literature 
presented by several researchers who investigated the physiology of the semicircular 
canals and the otolith organs, and also studied rotational and linear motion sensation. 
The development of the semicircular canals sensation model follows a previous 
presentation [13]. In addition, research on motion thresholds was surveyed in order to 
produce values to be used in the motion cueing algorithm development. 

The vestibular system is located in the inner ear and consists of the semicircular 
canals and otolith organs that sense angular and linear motion respectively. The location 
and orientation of the vestibular system in the head is shown in Figure 3.1. 



Figure 3.1. Location and Orientation of the Semicircular Canals. Reproduced with 
Permission from Purves, et al. [14]. 
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3.1. Semicircular Canals 


3.1.1. Physiological Description 

The semicircular canals consist of two sets of three elliptical cavities or canals 
that are filled with a fluid known as endolymph. The orientation of the canals in the head 
is shown in Figure 3.1. With the head in its normal erect position, the plane through the 
diameter of each horizontal canal is inclined about thirty degrees above an earth- 
horizontal plane. The posterior vertical canal lies in an almost vertical plane, forming a 
45-degree angle with the frontal plane of the head. The anterior canal is also at a 45- 
degree angle with the frontal plane, forming a right angle with the posterior canal. 

At one point on each canal, the canal cavity swells to form a bulbous expansion 
called the ampulla that contains the sensory ephithelium or crista. The crista contains 
bundles of sensory hair cells that extend into a gelatinous mass called the cupula as 
shown in Figure 3.2. The cupula bridges the width of the ampulla cavity, forming a seal 
through which endolymph cannot circulate. When the head turns in the plane of one of 
the canals, the inertia of the endolymph produces a force across the cupula, deflecting it 
in the opposite direction of head movement and causing a displacement of the hair 
bundles in each hair cell. Each hair cell has about 70 stereocilia and one kinocilium [15], 
with the stereocilia graded in length towards the kinocilium. Within one cupula, each 
kinocilium is on the same side as its stereocilia, forming a direction of polarization. 
When the cupula deflection is in the direction of the kinocilium, the hair cells will be 
maximally excited; whereas when the deflection is in the opposite direction the cells will 
be maximally inhibited. 
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Figure 3.2. Physiology of the Semicircular Canals Cupula. Reproduced with 
Permission from [14]. 

There are two types of sensory cells located in the cupula. The type I cells are 
contained in a nerve chalice and are innervated by fibers with a large diameter. The type 
II cells are cylindrical and are innervated by fibers with a small diameter. Both types of 
cells have a series of small hairs that penetrate into the cupula mass. 

3.1.2. Mathematical Model 

Zacharias [16] reported that Steinhausen first developed a linear second-order 
model of canal dynamics to explain the observed characteristics of vestibular-induced eye 
movements in fish (pike). This model was further refined by the “torsion-pendulum” 
model of Van Egmond, et al., [17], and was later developed from a systems approach by 
Mayne [18]. The transfer function for this overdamped system is 

hh , (3.D 

a(s) ( 1 + 7 1 5')(1 + t 2 s ) 

Further studies showed that the torsion-pendulum model does not completely 
represent rotational sensation. Young and Oman [19] formulated an adaptation operator 
and cascaded it with the torsion-pendulum model to resolve the conflicts between the 
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response predicted by the torsion-pendulum model and the perceptual responses 
measured in experiments. The addition of the adaptation operator resulted in the 
following transfer function: 


<t>c{s) _ K I" T a s ¥ 1 

a(s) scc [(1 + T a s) J|_(l + ^)(1 + t 2 s) 


(3.2) 


where the gain K scc noted by Zacharias [16] is proportional to ^ r 2 . 

Zacharias [16] reported several experiments suggesting an additional lead 
component. With the addition of this component, a model representing both the 
semicircular canal dynamics and the neural transduction dynamics was established: 


tc( s ) _ K \ v ir ( i+ ^) 

oc(s) SCC |_( 1 + V)J |_( 1 + V)(l + V)_’ 


Parameters for man are difficult to measure because direct measurement of the 


afferent response of the semicircular canals cannot be obtained. Therefore, most early 
experiments to determine the torsion-pendulum model parameters were based on 
subjective responses. Van Egmond [17] reported that the long time constant T t and short 
time constant r 2 had values of about 10 seconds and 0.1 seconds respectively. The values 
were based on the verbal response of humans subjected to various motion inputs in both a 
rotating chair and a torsion swing. Zacharias [16] noted that Meiry, measuring detection 
latency as a function of angular acceleration step size, obtained a 7- second long time 
constant for roll-axis rotation about the earth- vertical axis, and that Guedry, using a short 
period rotational stimulus consisting of an acceleration pulse doublet, and a response 
measure of apparent displacement, found values of 16 seconds for yaw-axis rotation and 
7 seconds pitch-axis rotation about the earth-vertical axis. Zacharias [16] then reported 
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that Malcolm and Melvill-Jones investigated the response to earth-vertical rotation about 
all three axes by using a velocity step as the stimulus, and measured the elapsed time to 
zero perceptual response. They obtained values of 6.1 seconds for the roll axis, 5.3 
seconds for pitch, and 10.2 seconds for yaw. 

Goldberg and Fernandez [20] determined average parameters for the semicircular 
canals of the squirrel monkey by direct measurement of the afferent nerves due to various 
angular acceleration inputs of different amplitudes and frequencies. Their transfer 
function related the afferent firing rate of the vestibular nerve to the angular acceleration 
input: 

AFR(f) [ 80s T (1 + 0.049,) 1 

ff(j) ' [(1 + 80 j)J[(1 + 5.7j)(1 + 0.003j)J 

The model parameters were estimated with the exception of the short time 
constant T,, which was determined analytically based on the physiology of the 
endolymph. Goldberg and Fernandez [20] noted that the short time constant r 2 is 
estimated to be 0.005 seconds for humans. 

It can be inferred that the long time constant r, measured by Van Egmond, et al., 
Meiry, Guedry and Malcolm and Melvill-Jones, as reported by Zacharias [16], does not 
actually represent the semicircular canal parameter in the model, but is an overall 
dynamics parameter representing the rotational sensation response to an angular velocity 
input. Zacharias [16] suggested that each axis of rotation has an equivalent “body axis” 
canal pair with a distinct time constant. The psychophysical results show each of the 
three axes having a distinct value for T x .. However, physiological results based on 
afferent responses by Goldberg and Fernandez show the same value for r, for the three 
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canal pairs. Zacharias [16] suggested the differences shown in the psychophysical results 
could occur at a central origin at the perceptual level. 

From subjective pilot measurements of angular acceleration thresholds on a 
moving base platform, Hosman and Van der Vaart [21] obtained the following 
semicircular canals transfer function, neglecting gain sensitivity and adaptation: 


■(*) = 


1 + 0.10975 


(1 + 5.924s)(l + 0.005s)' 


These results are based upon roll and pitch acceleration thresholds; yaw 
thresholds were not measured. The value for % agrees well with the value obtained by 
Goldberg and Fernandez [20]. The value obtained for t l compares to a value of 0.06 
seconds that Zacharias [16] reported that Benson and Ormsby obtained in experiments 
measuring nystagmus or involuntary eye movement due to motion. 

Zacharias [16] assumed that the angular velocity 6) from the semicircular canals 
that is sensed by human subjects is proportional to the cupula deflection <fc, and is 
expressed by the transfer function 


cb(s) _ Tj5 
co(s) (1+ r,5)(l + T 2 5)’ 


(3.6) 


where Zacharias [16] noted that the sensitivity gain is equal to the magnitude of the long 
time constant T h Goldberg and Fernandez [20] obtained gain sensitivity between the 
input stimulus and the afferent firing rate that was estimated at 3.44 spikes/sec per 
deg/sec 2 . Zacharias [16] noted that Ormsby suggested that the sensed angular velocity co 
is proportional to the afferent firing rate. While no one to date has experimentally 
obtained this parameter, Zacharias [16] reported that Curry, et al. provided an estimate of 
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the overall gain between perceived and input angular velocity based on angular 
acceleration thresholds. 

From this research, a transfer function that can best relate the sensed angular 
velocity to the acceleration stimulus is employed in the motion cueing algorithm 
development: 


®(' s ')_ 5 73 805 (1 + 0.06$) 
a(s) ~ ' (1 + 80$) (l + 5.73$)(l + 0.005$) ' 


(3.7) 


The frequency response of the transfer function given in Eq. (3.7) is shown in 
Figure 3.3. Both the torsion-pendulum model and the complete model with the lead and 
adaptation mechanisms included are shown. 


Semicircular Canals Transfer Function 
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Figure 3.3. Frequency Response of Semicircular Canals Transfer Function. 


The sensory function of the semicircular canals can be described by observing the 
frequency response of the torsion-pendulum model. In the range of normal head 
movement from 0.1 to 1.0 Hz [22], the gain response decreases by 20 dB/decade with the 
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phase close to minus 90 degrees. In this frequency range, the semicircular canals 
function as “integrating accelerometers” or angular velocity transducers. At very low 
frequencies less than 0.01 Hz, the phase approaches zero degrees, thus functioning as an 
accelerometer. At very high frequencies greater than 100 Hz, the phase approaches 
minus 180 degrees, thus functioning as an angular displacement transducer. The effects 
of adaptation and lead on rotational sensation are apparent; adaptation influences the 
afferent response at low frequencies below 0.01 Hz while the lead component influences 
high frequencies greater than 10 Hz. 

For implementation into the optimal cueing algorithm, angular velocity is 

employed as a stimulus, requiring the following transfer function: 

^( 5 )_ 5 73 8Q ^ (1 + 0-065) 

co(s) ‘ (1 + 80s)(l + 5.73s )(1 + 0.0055)' 

In addition, numerical stability problems may result when integrating the transfer 
function due to the small magnitude of the short time constant t 2 in the denominator 
relative to the simulation time step. Solely neglecting the short time constant would 
result in an unrealizable transfer function, but the lead time constant t l in the numerator 
could also be neglected since its order of magnitude is the same as the cueing algorithm 
time step. For numerical integration, the step size should be at least ten times smaller 
than the smallest time constant. The effect of both r 2 and T, is also well above the range 
of normal head movements. For these reasons a reduced-order transfer function can be 
utilized: 


g(f) 5 „ 8(fa 

a(s) ‘ (l + 80s)(l + 5.735) 


(3.9) 
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3.1.3. Physiological Interpretation 

Wu [13] presented a physiological interpretation of the behavior of the 
semicircular canals. He noted that the cupula-endolymph system first transforms an 
acceleration input to the head into a displacement of the cupula. This displacement then 
becomes an afferent response through a “mechano-neural” transduction system consisting 
of sensory hair cells and both efferent and afferent nerves. 

Many researchers have shown that an overdamped torsion-pendulum model could 
represent the cupula-endolymph system. Wu [13] reported that the remaining terms 
represent an adaptation-lead mechanism, noting the controversy over whether its origins 
lie in either the cupula-endolymph or the mechano-neural system. Wu reported that 
Goldberg and Fernandez [20] assumed that the origin of the adaptation mechanism might 
be centered on the physiology of the hair cells and/or the afferent neurons. Wu then 
noted that Goldberg and Fernandez suggested the lead mechanism may originate from 
sensory hair cells that are sensitive to both cupula displacement and velocity, which is 
reflected in the time constant r L . Wu [13] presented an interpretation by Schmid, et al. in 
which the lead mechanism is represented by efferent pathways that modify the feed- 
forward afferent dynamics by means of a negative feedback. Wu [13] demonstrated that 
this approach would justify the difference in order-of-magnitude of the adaptation time 
constant r A and the lead time constant r L . 

3.2. Otoliths 

3.2.1. Introduction 

The otolith organs are the elements of the vestibular system that provide linear 
motion sensation in humans and mammals. These organs are responsive to specific 
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force, responding to both linear acceleration and tilting of the head with respect to the 
gravity vector. However, the otoliths cannot discriminate between acceleration and tilt, 
requiring additional sensory information to resolve this ambiguity. There are two otolith 
organs, the utricle and saccule, in each inner ear. The utricle primarily senses motion in 
the horizontal plane, while the saccule primarily senses motion in the vertical plane. The 
otolith organs are inclined by about 20 to 30 degrees above the earth-horizontal plane as 
shown in Figure 3.1. 

3.2.2. Physiological Description 

The otolith organs consist of a two-layer structure known as the otolithic 
membrane that is attached to a base containing sensory cells. The otolithic membrane is 
composed of an upper layer, the otoconial layer, and a lower layer, the gelatinous layer. 
The endolymph fluid is in contact with the upper surface of the otoconial layer. The 
otoconial layer consists of calcium carbonate crystals embedded in a gelatinous material 
that rests on a less dense and extremely deformable gelatinous layer. This gelatinous 
layer is in turn attached to the sensory cell base known as the macula that is incorporated 
into the membranous tissue walls of the inner ear. The macula is rigidly attached to the 
skull and therefore moves with the head. 

There are two types of sensory cells located in the macula. The Type I cells are 
enclosed in a nerve chalice and are innervated by fibers with a large diameter. The Type 
II cells are cylindrical and are innervated by fibers with a small diameter. Fernandez and 
Goldberg [23] reported that cells in the outer (peripheral) otolith region are primarily 
Type II cells, and in the central (striolar) region cells are primarily Type I. Both types of 
cells have a series of small hairs that penetrate the lower portion of the gelatinous layer. 
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Each hair cell has about 70 stereocilia and one kinocilium [15], with the stereocilia 


graded in length toward the ki nocilium. 

The resulting displacement of the otolithic membrane due to forward linear 
acceleration is illustrated in Figure 3.4. The arrow in the figure shows the direction of the 
specific force acting upon the head. With a forward acceleration or backward tilting of 
the head, the denser otoliths tend to lag behind the macula, with the relative motion 
resulting in deformation of the gelatinous and otoconial layers in shear. When the shear 
deformation is in the direction of the kinocilium, the cell will be excited, whereas when 
the deformation is in the opposite direction, the cell will be inhibited. The directions of 
the maximum excitation and inhibition of a hair cell are defined by its polarization axis. 
In each macula, the central parting known as the striola separates oppositely polarized 
regions. For each position due to translational movement, some cells will be maximally 
excited, while others will be maximally inhibited. 


Figure 3.4. Displacement of the Otolithic Membrane due to Forward Acceleration. 
Reproduced with Permission from Purves, et al. [14]. 

Fernandez and Goldberg [24] identified two types of neurons that are 
characterized by their variance or regularity of discharge, hereafter referred to as regular 
and irregular units. From a sample population of units, they identified a ratio of regular 


Acceleration 
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to irregular units to be approximately three to one. The response of a neuron is the 


afferent firing rate (AFR). 

3.2.3. Mathematical Modeling 

Zacharias [16] reported that Meiry first investigated subjective responses to linear 
motion by using a cart to produce longitudinal sinusoidal motion. By measuring the 
subjective indication of direction, he obtained a transfer function relating perceived 
velocity v to stimulus velocity v: 

v(*) 


K T V 

OTO M * 5 


v{s) (v + l)(v + l)’ 


(3.10) 


where the long time constant T x and short time constant r 2 are 10 and 0.66 seconds 
respectively, and the gain K OTO is undetermined since amplitude measurements were not 
taken. Zacharias [16] then noted that Peters suggested the subjective response measured 
by Meiry was perceived acceleration and not perceived velocity, since in response to an 
acceleration step, the model predicted a perceived response that decays to zero with a 
time constant of 10 seconds. 

Young and Meiry [25] noted that the model proposed by Meiry correctly 
predicted the phase of perceived velocity for lateral oscillation and time to detect motion 
under constant acceleration, but failed to predict the otoliths’ response to sustained tilt 
angle as indicated by behavioral and physiological data. They noted that the model 
agreed with dynamic counter-rolling data (of the eye) at high frequencies, but 
experimental counter-rolling at zero frequency showed a static component of otolith 
output with no phase lag (the model assumed no static output and at zero frequency 
approached 90 degrees of lead). They proposed the following revised model of specific 
force sensation: 
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(3.11) 


f{s) 1.5 (5 + 0.076) 


f{s) (s + 0.19)(s + 1.5)’ 


which, when rearranged in terms of the time constants, yields 


f(s) 0.4(13.25 + 1) 


f(s) {5.33s + l)(0.66s + l)' 


(3.12) 


With a smaller long time constant (5.33 seconds) and an additional lead term, they 
modeled both perceived tilt and acceleration in response to acceleration input. They 
noted that the model acts as a velocity transducer over the frequency range of 0.19 to 1.5 
rad/s, with the transfer function from specific force to perceived tilt or lateral acceleration 
having a static sensitivity of 0.4. This model presumes the equivalence of linear 
acceleration sensation with that of tilt. 

Zacharias [16] noted that a lumped parameter model of otolith motion could be 
used to represent the two lag time constants, similar to the torsion-pendulum model for 
the semicircular canals. Ormsby [26] first developed this model, and Grant, et al. [27-30] 
later refined the model as part of their theoretical analysis of the otolithic membrane. 
Grant and Best [30] obtained the following transfer function for the model: 


*(s) 

m 


I-* 


T \ T 2 


aJ(!+ V)( 1 + V) 


(3.13) 


where x is the relative displacement of the otoconial layer with respect to the head, p e is 
the density of the endolymph, p 0 is the density of the otoconial membrane, with p„ > p e . 
For the otoliths, we again have an overdamped system with x, » x 2 . 

In determining the value of the short time constant r 2 , Grant and Best [30] first 
examined the maximum displacement of the otoconial layer in response to a step change 
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in linear velocity. The acceleration for a linear velocity step U is a x = —U S(t ) , with g x = 
0, where 8(t) is the unit impulse function. The transient response to Eq. (3.13) is then 

( p \ -z -z 

x(t) = U 1-^ T 2 (e A -e /Tl ). (3.14) 

v Ay 

By assuming that the short exponential term in Eq. (3.14) has reached zero and 
the long exponential term remains close to unity, the maximum displacement of the 
otoconial layer v max can be approximated as 



The theoretical continuum mechanics analysis performed by Grant and Best [29] 
first indicated that this short time constant t 2 is 0.002 seconds or less. They later 
demonstrate that this value turns out to be too large when reasonable values of the 
maximum otolith displacement are considered. For p„ = 2.0 and U = 25 cm/sec (a 
reasonable value for normal head velocity), Eq. (3.15) becomes v max = 12.5 r 2 . For r 2 = 
0.002 sec, the maximum displacement of the otolithic membrane resulted in v max = 250 
pm. It was assumed for shear deformation the maximum displacement should not exceed 
the thickness of the otoconial layer (25 pm), indicating the short time constant should be 
one order of magnitude smaller, r 2 = 0.0002 sec. This indicated that more damping was 
needed in the lumped parameter model. Grant and Best [30] showed that additional 
damping could be introduced by the inclusion of a viscoelastic gelatinous layer in the 
continuum mechanics model. 

Ormsby [26] neglected the short time constant t 2 in Eq. (3.15), and after 
rearranging terms, approximated the otolith mechanical dynamics by 
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(3.16) 


x(s) _ A 


and then proposed a model for the response of the otolith afferent dynamics: 

AFR(s) _ Bs + (B + C)A 


m 


s + A 


(3.17) 


where AFR is the change in afferent firing rate from the resting discharge, and the 
constants A, B. and C are undetermined. This model assumes that higher centers process 

the afferent response optimally to estimate the perceived specific force / as shown in 
Figure 3.5. 


/w l 

Bs + (B + C)A 

AFR(sY 

H( S ) 

/(*) 

► 

s + A 


► 


Combined Processing by 

Mechanical and Higher Centers 

Afferent Otolith 
Dynamics 

Figure 3.5. Model of Otolith Specific force Sensation. Adapted from Ormsby [26]. 


The steady-state optimal processor H(s) is then determined by solving the 
associated Wiener-Hopf equation [31], yielding a solution of the form 

5 + A 


H(s) = M 


(s + F)(s + Gy 


(3.18) 


where F, G, and M are nonlinear functions of the independent variables A, B, and C in 
Eq. (3.17). H(s) is then cascaded with the otolith afferent dynamics to estimate the 
perceptual dynamics associated with the otoliths: 
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(3.19) 


l 

f 



which is equivalent to the Young-Meiry [25] model given in Eq. (3.11). Ormsby [26] 
noted that Fernandez, Goldberg, and Abend found an average steady-state change in 
afferent firing rate from the utricle due to a 1-g step to be 45 impulses per second (ips), 
resulting in the condition that B + C = 45. Setting Eq. (3.19) equal to Eq. (3.11) and 
including this constraint results in the following model for the afferent response: 


af/?co _ 9Q j + o-1 

m .v+o.2’ 


(3.20) 


This transfer function, when rearranged in terms of its time constants, becomes 


AFR(s) _ 105 + 1 

f ( 5 ) “ 55 + 1' 


Ormsby [26] noted the following about this model: 

The approach taken here can yield a model which accounts reasonably 
well for the available subjective data, the known physiological structure of 
the sensor and makes reasonable predictions concerning the afferent 
processes and the associated central processing. 

Fernandez and Goldberg [23] studied the discharge of peripheral otolith neurons 
in response to sinusoidal force variations in the squirrel monkey. Both regularly and 
irregularly discharging neurons were measured. The gain curves for the regular units 
were flat, with a small phase lead at low frequencies and a larger phase lag at higher 
frequencies. The irregular units showed a larger gain enhancement and phase lead at 
high frequencies, which could not be represented by a first-order lead operator. They 
noted on average, there is an increase by a factor of eighteen in gain enhancement in 
irregular units but only an increase of a factor of two for regular units. 
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The frequency responses of regular and irregular units resulted in a transfer 
function of the form 


1 + k A r A s 1 + k v ( T v s) k 


1 + Z A S 




In Eq. (3.22), the term H v is a velocity- sensitive operator with a fractional exponent (k v < 
1) and provides most of the gain enhancement and phase lead found in the irregular units. 
The value of k v reflects the effectiveness of the lead operator and is closely related to the 
slope of the gain curve. The term H A is an adaptation operator that contributes to low 
frequency phase leads and increases of gain from static or zero frequency to 0.006 Hz. 
The term H M is a first-order lag operator that Fernandez and Goldberg [23] noted might 
reflect the mechanics of otolith motion. This lag term accounts for the high frequency 
phase lags observed in regular units and for high frequency phase leads in irregular units 
being smaller than would be predicted solely by a fractional lead operator. The term K 0T0 
defines the static sensitivity in terms of afferent firing rate per unit of acceleration, i.e., in 
units of impulses per second per g (ips/g). 

Fernandez and Goldberg [23] estimated parameters for the transfer function, 
obtaining nearly equal results for various values of T v . The median parameters for both 
regular and irregular units for T v = 40 seconds are given in Table 3.1. 


Table 3.1. Median Parameters for Regular and Irregular Units. 



k v 

k A 

Ta 


Koto 

Regular 

0.188 

1.12 

69 sec 

16 msec 

25.6 ips / g 

Irregular 

0.440 

1.90 

101 sec 

9 msec 

20.5 ips / g 


Because of the fractional exponent in the transfer function of Eq. (3.22), an 
elementary solution to its response cannot be readily obtained. However, an approximate 
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solution to the response can be derived through the application of fractional calculus [32]. 
The derivation to obtain a response to a transfer function with fractional exponents is 
given in Appendix A. Applying this derivation with the regular unit parameters given in 
Table 3.1 results in the response to a unit step: 

x(t) = 25.601 - 28.673A 62 ' 5 ' + 3.073^ 0014493 ' 


10.786E, (-0.188,-62.5) + 1.156E, (-0.188,-0.014493) (3.23) 

f 


+ 9.629 


r'O'+i)' 


where T is the gamma function. 

Eq. (3.23) is an infinite series. For v equal to zero, Eq. (3.23) will reduce to the 
Taylor series expansion of the exponential function. When v is not equal to zero, 
E t (v, a) is a transcendental function that can only be approximated. A recursion 


formula was derived, where the solution to the function E t (v, a) is given as 


E,(v,a) 
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at 
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r(v + i) r(v + 2) r(v + 2) 


I‘ 


-au v+\ i 

e u au. 


(3.24) 


Similarly, the unit step response for the irregular unit parameters is derived: 


x(t) = 20.308 - 35.588A 111 ' 1111 ' + 18.280 A 0009901 ' 

- 86.063E, (-0.44, -111.1111) + 40.769E, (-0.44, -0.009901) 


(3.25) 


+ 45.294 


r ( |/+1 ) 


Hosman [33] noted that the fractional exponent models are not easy to implement 
in motion cueing algorithms due to the fractional exponent in the lead term. He reported 
a simplified model of the same form developed by both Ormsby [26] and Grant and Best 
[29]: 
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(3.26) 


AFR(s) _ 33 3 (* + 1) 

f(s) ' (0.5s + l)(0.016s + l)' 

Note that the gain terms for the Femandez-Goldberg model [23] from Table 3.1 
are about one half that of the gain value used by Ormsby [26] to develop his model. 
Hosman [33] proposed a gain term of less magnitude than that used by Ormsby [26] that 
may provide an improved approximation to the Fernandez-Goldberg responses. Due to 
the adaptation mechanism in the Femandez-Goldberg model [23], these gains will require 
a long duration step input to be realized in steady state. Hosman [33] chose the short 
time constant r 2 to be equal to the otolith mechanics time constant t m reported by 
Fernandez and Goldberg [23] for the regular units. No basis, however, was given for the 
values selected for the long time constant % and the lead time constant t l , which are one 
order of magnitude less than those resulting from the model developed by Ormsby [26]. 

By using the long and lead time constants reported by Ormsby [26] in Eq. (3.21), 
and selecting the short time constant and gain reported by Hosman [33] in Eq. (3.26), the 
following transfer function results for the afferent otolith dynamics: 

^1 = 333 < 10 * + 1 ) (3 27) 

/'III ■ (5j + 1)(0.016j + 1)' 

The response to a step input of 1-g magnitude (9.81 m/s“) will now be examined 
for the Fernandez-Goldberg model [23] with both the regular and irregular unit 
parameters. Figure 3.6 compares the step responses to the response for the proposed 
afferent dynamics model given in Eq. (3.27) for both 1 -second and 30-second durations. 
Note that the onset for the proposed model is faster than the regular unit, but slower than 
the irregular unit. There is no large overshoot as observed with the irregular unit 
response. The steady-state response for the proposed model is less than the irregular unit 
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response but greater than the regular unit response, and approaches the regular unit 
response for the given time duration. Both the regular and irregular unit response will 
slowly approach their respective gain values, and beyond about 80 seconds the irregular 
unit response will decrease below that of the proposed model. The model more closely 
represents the population-dominant regular units with a faster onset and higher magnitude 
steady state effects that occur in the less prevalent irregular units. 


Comparison of Otolith Models Response to Step Input 



0 0.5 1 0 10 20 30 

Time (sec) Tme (sec) 

Figure 3.6. Comparison of Otolith Models Response to a 1-g Step Input. 


From this research, a transfer function that can best relate the sensed response to 
the specific force stimulus is proposed: 


/ , (Mill 

/ oro (v + 1X^+1)’ 


(3.28) 


with K oto = 0.4, r, = 5 sec, t 2 = 0.016 sec, and r L = 10 sec. For implementation into the 
motion cueing algorithms, Eq. (3.28) can be rewritten as 


/ , (j + A,) 

/ oro (s + 5 0 )(s + A)’ 


(3.29) 
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where A 0 = 1 /t l , B„ = 1 / r, . B, = l/r 2 , and K' oro = K oto tj 2 /t l . 


The frequency response of the proposed transfer function in Eq. (3.28) compared 
to the Young-Meiry model [25] of Eq. (3.11) is shown in Figure 3.7. Note that the gain 
and phase lag for the Young-Meiry model [25] occur at a much lower frequency as 
compared to the proposed model. This is due to the magnitude of the short time constant 
t 2 for the Young-Meiry model [25] being an order of magnitude larger than the value 
used in the proposed model that was obtained by Fernandez and Goldberg [23]. In the 
range of normal head movements from 0.1 to 1.0 Hz [22], the gain for the proposed 
model remains constant, with the phase close to zero degrees. In this frequency range, 
the otolith functions as a specific force transducer. 


Frequency Response of Otolith Models 



Figure 3.7. Frequency Response of Proposed and Young-Meiry Otolith Models. 
3.2.4. Physiological Interpretation 

Modern theories of the operation of the otolith receptors are based on the 
assumption that the afferent responses are generated by the deflection of hairs in the 
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sensory cells as a result of the otolith displacement. A specific force input, in the form of 
either linear acceleration or tilt, is transformed by the otolith-endolymph system into a 
displacement of the otolith. This displacement is further transformed into an afferent 
response by the mechano-neural transduction system consisting of sensory hair cells and 
both efferent and afferent nerves. 

Many researchers have shown that the otolith-endolymph system could be 
represented by an overdamped mass-spring-damper system. Grant and Best [29] reported 
that the magnitude of the long time constant % is considered correct by most investigators 
because the overall system (otolith organ, neural transmission, central nervous system 
processing, and involuntary eye motion) could easily follow such a slow time constant. 
The value Grant and Best [29] obtained for the short time constant r 2 is a two order-of- 
magnitude decrease in time constant as compared to the value obtained from the ocular 
torsion responses measured by Young and Meiry [25]. This value is also a one order-of- 
magnitude decrease as compared to the value of r M that Fernandez and Goldberg [23] 
attribute to the afferent dynamics. The fast dynamic response of the otolith will decrease 
to the slower ocular torsion response due to losses in neural transmission and central 
nervous system processing. 

Young and Meiry [25] first noted that the origin of the lead term could be 
neurological, either in central processing of the otolith displacement signals or through 
the presence of two types of hair cells in the macula. One type of hair cell would respond 
to displacement and the other would respond to the rate of change of otolith 
displacement. These hair cells could produce the lead term if they were of the slowly 
adapting type postulated by several researchers. Fernandez and Goldberg [23] later show 


44 



that the degree of sensitivity to the otolith velocity is a function of the fractional 
derivative in the lead operator, i.e., irregular units are more velocity- sensitive than regular 
units. They noted that this difference in sensitivity might be due to discrepancies that are 
more noticeable in irregular units. 

Fernandez and Goldberg [23] suggested that the difference between the expected 
otolith displacement and afferent firing rate for both regular and irregular units may be 
attributed to the mechanical linkages between the sensory hair bundles and the gelatinous 
layer. They reported that these sensory hair bundles are not rigidly embedded in the 
membrane, but are enclosed in a fluid-filled meshwork between the membrane and the 
sensory epithelium. Motion could be transferred to the hairs either by directly contacting 
the meshwork or indirectly by viscous coupling with the fluid. They also noted that the 
irregular units correspond to thick afferents that stimulate the Type I hair cells in the 
striola. Grant and Best [30] also suggested that the nonlinear stiffness of the gelatinous 
layer could also contribute to these differences as well. 

3.3. Motion Thresholds 

Zacharias [16] reported that Clark reviewed twenty-five earlier studies that 
attempted to define an absolute threshold for angular acceleration. Clark noted the wide 
range in rotational devices, stimuli waveforms, psychophysical methods, and threshold 
definitions employed by various researchers. The threshold measurements reported 
showed a two order-of-magnitude range for yaw-axis earth-vertical rotation. Zacharias 
also noted that of the twenty-five studies that Clark reviewed, only one study was 
reported for pitch and one for roll rotation. 


45 



Clark and Stewart [34] later conducted a study with angular acceleration step 
input stimuli. Using fifty-three test subjects, they found a mean threshold for the 
perception of the yaw-axis earth-vertical rotation to be 0.41 deg/sec“. In a separate 
study, Clark and Stewart [35] studied thresholds about all three earth-vertical axes for 
eighteen subjects. The mean threshold for yaw and roll was about the same (0.41 
deg/sec 2 ), and was found to be larger for pitch (0.67 deg/sec 2 ). 

Zacharias [16] reported that Mulder first recognized that the product of 
acceleration magnitude with detection or latency time (Mulder product) is approximately 
a constant, thus suggesting the existence of an angular velocity threshold. Zacharias [16] 
noted that Van Egmond demonstrated that the Mulder product could be derived from the 
torsion-pendulum model, resulting in an estimated value of about 2 deg/sec. Meiry [36] 
employed a step-response technique to measure latency time of subjects in a motion 
simulator. For yaw-axis earth-vertical rotation, Meiry [36] obtained a value of 2.6 
deg/sec. For roll-axis earth- vertical rotation, a value of 3.0 deg/sec was obtained. 

Zacharias [16] then demonstrated that given an infinite detection time 
corresponding to acceleration below an absolute threshold, the velocity threshold is equal 
to the acceleration threshold multiplied by the semicircular canals long time constant. He 
then showed that multiplying the measured accelerations obtained by Clark and Stewart 
[35] by the long time constants for the corresponding axis obtained by Melvill-Jones, et 
al. resulted in estimated angular velocity thresholds for each earth-vertical body axis. 
Zacharias noted that these values were in general agreement with those obtained by 
Meiry [36]. 
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Benson, et al. [37] reported that threshold measurements were never made with an 
angular motion stimulus having a trajectory similar to that of a natural head movement. 
To address this, they performed threshold experiments on a turntable driven by a 
precision torque motor that generated rotational stimuli about an earth-vertical axis, 
carrying either a seat with the subject seated upright, or a stretcher that allowed the 
subject to lay either supine or on their right side. The following results in Table 3.2 were 
obtained for all three earth-vertical axes, and are compared to the thresholds reported by 
Reid and Nahon [9] that include the roll and yaw thresholds reported by Meiry [36] and 
the pitch threshold Zacharias [16] estimated from the acceleration thresholds obtained by 
Clark and Stewart [35]. 


Table 3.2. Comparison of Body Axis Angular Velocity Thresholds in deg/sec. 


Reference 

Roll 

Pitch 

Yaw 

Reid and Nahon 

3.0 

3.6 

2.6 

Benson, et al. 

2.0 

2.0 

1.6 


Zacharias [16] reported a review of linear acceleration threshold studies by Peters, 
in which Peters noted a one order-of-magnitude range in measured threshold (0.002 to 
0.02 g). Possible contributions to this variation included the variability between subjects, 
the type of stimulus used (e.g., sinusoidal vs. step), the definition of threshold, and the 
head axis orientation with respect to the stimulus. Zacharias [16] noted that only one of 
the reviewed studies used a linear acceleration stimulus in the earth- vertical direction, in 
which Mach obtained an acceleration threshold of 0.012 g. Subsequent vertical motion 
threshold measurements reported by Zacharias [16] show almost a one order-of- 
magnitude difference, from 0.0085 g to 0.005 g. 
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Zacharias [16] also reported threshold measurements with a linear acceleration 
stimulus in the horizontal plane with subjects lying supine, noting that similar thresholds 
to vertical motion are expected since this geometry implies an alignment of the stimulus 
acceleration vector with the vertical body axis. Zacharias [16] noted that these results are 
confirmed by an estimate of 0.0 lg obtained by Meiry [36] (using a linear motion cart and 
acceleration step inputs). However, Zacharias [16] then noted that the most common test 
protocol for threshold measurements of horizontal stimuli has been with subjects seated 
upright. Meiry [36] noted that with the utricle inclined at about 30 degrees above the 
horizontal head plane, horizontal thresholds might be expected to be lower than vertical 
thresholds by a factor of cos 30/sin 30, or about a factor of 1.7. 

Benson, et al. [38] performed threshold experiments with a test stimulus 
consisting of a single discrete movement having an acceleration trajectory that 
approximated a sine wave. This stimulus is similar to the trajectory used by Benson, et 
al. [37] in determining rotational thresholds. Motion stimuli were generated by a 
horizontal linear oscillator guided by externally pressurized aerostatic bearings, 
supporting a seat assembly that could be adjusted so that the stimuli axis was parallel to 
the axis of motion of the carriage. Benson, et al. [38] obtained thresholds for the x-, y-, 
and z-axes that are noted in Table 3.3, and are compared to thresholds reported by Reid 
and Nahon [9]. Reid and Nahon [9] reported acceleration threshold values that were 
based on the studies reported by Zacharias [16]. The z-axis threshold is about the same 
as that Zacharias [16] reported was obtained from the Mach study reported by Peters. 
The x- and y-axis acceleration is about a factor of 1.7 less than that noted for the z-axis, 
which is consistent with the observation noted by Meiry [36]. Benson, et al. [38] noted 
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that these thresholds are within the ranges reported in other studies; in particular the 
results obtained by Meiry [36], who employed a step acceleration stimulus. The 
significantly higher threshold for the z-axis was noted by Benson, et al. to be consistent 
with the findings of studies employing sustained oscillatory or step acceleration stimuli. 


Table 3.3. Body Axis Linear Acceleration Thresholds in m/sec 2 . 


Reference 

X-Axis 

Y-Axis 

Z-Axis 

Benson, et al. 

0.0625 

0.0569 

0.154 

Reid and Nahon 

0.17 

0.17 

0.28 


The linear and angular motion thresholds presented in Tables 3.2 and 3.3 
respectively are effective, or “indifference” thresholds that are more appropriate for the 
pilot-vehicle environment than absolute thresholds that result from the detection of a 
single task in an ideal laboratory environment. Zacharias [16] noted that higher 
indifference thresholds during active tracking are justified because of less attention given 
to motion cues due to workload. Gundry [2] showed an increase of 40% when the 
subject is loaded with an arithmetic task. Hosman and Van der Vaart [21] observed a 
similar increase in roll and pitch thresholds when their subjects were loaded with either a 
control task or an auditory discrimination task. 

Zacharias [16] noted that a latency dependence on angular acceleration is 
observed, and a velocity threshold model similar to angular motion can be proposed. 
Meiry [36] measured detection latencies as a function of linear acceleration step size. 
This model assumed a velocity threshold exists such that acceleration thresholds required 
T seconds to be detected. In response to a velocity ramp input, the model predicted a 
perceived velocity. For subjects seated upright, the model resulted in a linear velocity 
threshold of 0.02 g-sec or about 0.2 m/sec for longitudinal motion. Zacharias [16] noted 
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that Melvill- Jones and Young used a similar analysis of detection time and acceleration 
from both their own experiments and those of Meiry [36]. Based on Meiry’s data, they 
proposed a value of 0.22 m/sec for both longitudinal and vertical motion. 

The angular velocity and linear acceleration thresholds given in Tables 3.2 and 
3.3 are used in the development of the linear optimal algorithm discussed in Chapter 4. 
The linear velocity thresholds mentioned in the last paragraph are incorporated in the 
integrated human perception model discussed in Chapter 5 and implemented in the 
nonlinear motion cueing algorithm developed in Chapter 6. The integrated perception 
model also includes the angular velocity thresholds of Table 3.2. 
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4. Linear Optimal Motion Cueing Algorithm 

4.1. Problem Description 

In developing a set of linear washout filters, the problem is to determine a transfer 
function matrix W(s) that relates the desired simulator motion input to the aircraft input 
such that a cost function constraining the pilot sensation error (between simulator and 
aircraft) is minimized. The structure of this problem is illustrated in Figure 4.1. 


Aircraft 
States u A 



Simulator Pilot 


Figure 4.1. Linear Optimal Algorithm Problem Structure. 


A mathematical model of the human vestibular system is used in the filter 
development. The optimal algorithm generates the desired transfer functions W(s) by an 
off-line program, which are then implemented on-line. W(s) will relate the simulator 
motion states to the aircraft states by u s = W(s) x u A . The simulator states u s are then 
used to generate the desired motion platform commands. 

In the original development, the washout filters were applied in the pilot head 
reference frame. Reid and Nahon [9] noted that this frame selection was chosen to 
eliminate sensation cross-couplings that made the development of W(s) more 
complicated. Wu [13] demonstrated that this location of the center of rotation at the 
pilot’s head resulted in excessively large actuator extensions in some input cases. He 
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suggested that the optimal algorithm be reformulated in the simulator reference frame 
with the center of rotation located at the centroid of the simulator motion-base. 

The question has arisen as to what aircraft and simulator control inputs are the 
most appropriate for the optimal algorithm. The previous developments suggested a 
control input for either the longitudinal or lateral mode with linear acceleration and 
angular displacement as control inputs. Wu [13] developed an approach using linear 
acceleration and angular acceleration for the longitudinal mode. This approach showed 
advantages in controlling additional modes that were not available in the original 
development. In addition, since the semicircular canals behave as a transducer for 
angular velocity input in the range of normal head movements, an approach using angular 
velocity as an input is desired. In this research an optimal algorithm based on simulated 
aircraft angular velocity inputs is developed. 

4.2. Algorithm Development 

4.2.1. Longitudinal Mode 

The algorithm development with angular velocity input for the longitudinal 
(pitch/surge) mode is given below. The input u is formulated as 


u = 

~9~ 

— 

u i 


a x _ 


_u 2 _ 


where 0 is angular velocity, and a x is the translational acceleration, with each term 
respectively set equal to m, and u 2 . 

The sensed rotational motion 9 is related to «, by the semicircular canals model 
of Eq. (3.8): 
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(4.2) 


GsccW ^ 

(1 + v)(l + v)(l + ^) 15 

where the semicircular canals time constants T u r 2 , T„, and t l are given in Eq. (3.8), and 
G scc is the angular velocity threshold that scales the response to threshold units. Eq. (4.2) 
can be rewritten as 



y \ s 3 + t 3 s 2 
s + 71,5" + ZJs + T q 


(4.3) 


where 
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1 o 


T _ s + r i + T - 
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. t,t, + r, (t, + r,) 

, T 2 = , T 3 = G scc / f 2 , and T 4 = G scc t l / f 2 , 

W2 


and can be defined in state space notation as 


x scc ^scc x scc ®scc u 


Q ^SCC X SCC ®SCC U ’ 


(4.4) 


where in observer canonical form, 
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C scc=[ 1 0 0], andD scc =[T 4 O], 


The sensed specific force f x is related to the stimulus specific force f x by the 
otolith model of Eq. (3.29): 


/,(■*) =g OT£ X 


(* + A>) 


OTO 


(s+ B 0 )(s + B { ) 


fx( S )’ 


(4.5) 


where A 0 , B 0 , B,, and K' oro are computed parameters of the otolith model as given in Eq. 
(3.29), and G 0T0 is the linear acceleration threshold that scales the response to threshold 


units. 
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For the center of rotation at the centroid of the motion platform, the specific force 


fx = a x+ ~ R s °- 


where R s . is the radius from the motion platform centroid to the pilot’s head. In terms of 
a, and u 2 , Eq. (4.6) is transformed into the Laplace domain, where 


( 1 ^ 

f x (s) = u 2 (s)+ g R Sz S Mj(s). 
V S 


Substituting (4.7) into (4.5) results in 

~/\ , ( s + A> ) f l 3 

fx( s ) - G 0 to R oto 7 TT7 |TT u 2 8 R Sz s u i 

(s + B 0 )(s + B^l y s J J 

, (s + Aq) 17 1 3 

= G oto K oto - , g R s - S 1 u 

(s + B oA s + 5 i)Lv s ; J 

= r K r \ ~R S / - R Sz A/ + 8 S + 8A) i S + A o) 1 

oro oro sis + B^s + B,) (s + B 0 )(s + Bl ) • 


Note that in Eq. (4.8) the system equation becomes realizable with the inclusion 
of the otolith break frequency 5, , which was neglected by both Reid and Nahon [9] and 
Wu [13] in their respective optimal algorithm formulations. Rearranging and taking 
derivatives results in the differential equation 

7+(s»+b,)A+w,= 

r . (4.9) 

Goto R oto R Sz ( R o "F ~ A) ) A "F ( 8 "F R sz R o R i )^i "F SAi y^\dt F tk + A () u 2 , 

which can be rewritten as 


/ v + af x + bf x = cu l + du { + e Jwjdf + fu 2 + gu 2 . 


(4.10) 


and can then be defined in state space notation as 
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v — A v -L R II 
^OTO ^OTO^OTO OTO 


/, = C c 


I D OTO U, 


where x OTO are the otoliths states, and 
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[10 0 1 0], D om =[-G OTO K' OTO r Sz 0]. 


(4.11) 


The representations in Eqs. (4.4) and (4.11) can be combined to form a single 
representation for the human vestibular model: 


x v = A v x v + B v u 

y 1 Cy Xy + Dyll, 


(4.12) 


where x v and y v are, respectively, the combined states and sensed responses, and A v , 


By, C v , and D v represent the vestibular models as one set of state equations: 
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It is assumed that the same sensation model can be applied to both the pilot in the 
aircraft and the pilot in the simulator as shown in Figure 4.1. We then define the 
vestibular state error x e = x s - x A (where x s and x A are the respective vestibular states for 
the simulator and aircraft), and the pilot sensation error e, resulting in 

x e =A v x e +B v u s -B v ii A 

C CyX e + Dytlg Dyll^ , 

where u s and u A represent the simulator and aircraft inputs as given in Eq. (4.1). 
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In order to constrain the simulator motion, it is necessary for the control algorithm 
to explicitly access states such as the linear velocity and displacement of the simulator 
that will appear in the cost function. For this purpose, additional terms are included in 
the state equations: 


X d =A d X d +B d U S’ 

where x d represents the additional simulator states: 


x d 





(4.14) 


and is related to the simulator input u s by 
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The aircraft input u A consists of filtered white noise, and can be expressed as 


x n = A „ X „ +B„w 


U A = X „> 


(4.15) 


where x„ are the filtered white noise states, w represents white noise, with A n and B n 


given as 


A 


n 


~~Yi 

0 " 

, B n = 


0 

~Yi_ 


7i_ 


where y t and y 2 are the first-order filter break frequencies for each degree-of-freedom. 
The state equations given in Eqs. (4.13), (4.14), and (4.15) can be combined to form the 
desired system equation 


x = Ax + Bu s + Hw 
y = [e Xj] 1 = Cx + Du s , 


(4.16) 
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where y is the desired output, and x = [x e x d x n ] T represents the combined states. 
The combined system matrices A, B, C, D, and H are then given by 



A cost function J is then defined as 

J = i?j (e T Qe + x d R d x d + UsRu s )dtj, (4.17) 

where E{ } is the mathematical mean of statistical variable, Q and R d are positive semi- 
definite matrices, and R is a positive definite matrix. Eq. (4.17) implies that three 
variables are to be constrained in the cost function: the sensation error e along with the 
additional terms x d and u s which together define the linear and angular motion of the 
platform. The cost function constrains both the sensation error and the platform motion. 

The system equation and cost function can be transformed to the standard optimal 
control form as shown in Kawkemaak and Sivan [39] and noted in Reid and Nahon [9] 
by the following equations: 

x = Ax + Bu' + Hw 

i' = £j J l (VR'x + u' 1 R 2 u')<r//j, ^ ^ 

where 

a' = a - br 2 'r' 2 , iT = u s + r^r'x, r; = r 2 - R 12 R 2 R[ 2 , 

R x = C t GC, R 12 = C t GD, R 2 = R + D'GD, G = diag[ Q, R d ] 

The cost function of Eq. (4.18) is minimized when 

u' = -R 2 1 B t Px, (4.19) 

where P is the solution of the algebraic Riccati equation 
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(4.20) 


R; - PBR^B'P + A ,t P + PA' = 0. 
Substituting Eq. (4.19) into Eq. (4.18) and solving for u s , 


u s =-[r- 1 (B t P + R 1 t 2 )]x, 


(4.21) 


and defining a matrix K, where u s = Kx , results in K = R ; 1 ( B 1 P + R J 2 j . 
K can be partitioned corresponding to the partition of x in Eq. (4.16): 


u s — [Kj K 2 K 3 ] x d 


(4.22) 


Noting that x n = u A , remove the states corresponding to the x n partition from Eq. (4.22): 


A v 0 -B v 
0 A d 0 


x d + R «s 


(4.23) 


and substituting Eq. (4.22) into Eq. (4.23) results in 


X e _ Ay-ByKj -ByK 2 X e ^ "By^ + Kj) ^ 

X d "^d^l ^d "®d^2__ x d_ _ ^d^3 


(4.24) 


After observing the state space form of Eqs. (4.24) and (4.13), the following equations 
are obtained in the Laplace domain: 


u s (s) = W(s)xu A (s), 


(4.25) 


where 


W(s) = [K, K.] 


Si - Ay + ByKj ByKj 


B d Kj sI-A d +B d K 2 J I B d K 3 


B V (I + K,) 


W(s) is a matrix of the optimized open-loop transfer functions linking the 


simulator inputs u s to the aircraft inputs u A . The block diagram for the on-line optimal 
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algorithm implementation is shown in Figure 4.2. Similar to the NASA adaptive 
algorithm [7], there are separate filtering channels for the translational and rotational 
degrees of freedom with the cross-feed path providing the tilt coordination cues. 

The aircraft acceleration input vector is first transformed from the aircraft body 
frame Fr A to the inertial frame Fq using Eq. (2.1). Nonlinear scaling in combination with 
limiting as described in Section 2.5 is then applied to scale the aircraft inputs. The scaled 
inertial acceleration a A is then filtered through the translational filter W 22 to produce a 
simulator translational acceleration command S, . This acceleration is integrated twice to 
produce the simulator translational position command Sj . 

The aircraft angular velocity input io A is transformed to the Euler angular rate 
vector P A using Eq. (2.2), and is limited and scaled similar to the translational channel. 
This input is then passed through the rotational filter W n to produce the vector P SR . The 
tilt coordination rate P ST is formed from the acceleration a A being passed through the tilt 
coordination filter W 12 . The summation of P ST and P SR yields P s , which is then 
integrated to generate P s , the simulator angular position command. 

The simulator translational position S 1 and the angular position P s are then used 

to transform the simulator motion from degree-of-freedom space to actuator space as 
given in Eqs. (2.6) and (2.7), generating the actuator commands required to achieve the 
desired platform motion. 
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Figure 4.2. Optimal Algorithm Implementation for Longitudinal Mode. 

The desired motion cueing filter matrix W(s) is computed using a set of 
MATLAB™ scripts. The weighting matrices Q, R, and R,i given in the cost function of 
Eq. (4.17) are selected and adjusted to produce the desired motion platform commands. 
From these weights and the vestibular models the standard optimal control matrices of 
Eq. (4.18) are computed. The algebraic Riccati equation of Eq. (4.20) is solved with a 
MATLAB™ function that uses a generalized eigenproblem formulation with a Newton- 
type refinement presented by Arnold and Laub [40]. The solution for W(s) is then 
computed. Common poles and zeros are cancelled in each transfer function, yielding a 
set of seventh-order filters for the longitudinal mode. These filters are then used in a 
SIMULINK™ model that generates the linear acceleration and angular velocity 
responses. If the solution to W(s) is unsatisfactory, this procedure is repeated by 
adjusting the elements of the weighting matrices Q, R, and Rd until the desired results are 
obtained. The procedure for computing the solution to W(s) is illustrated in Figure 4.3. 
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Compute System Matrices 
A, B, C, D from Vestibular, 
Platform, and White Noise 


Parameters 



Adjust Weighting Matrices 
Defined in Cost Function: 

IT R d > Q 



Transform System Equation 
and Cost Function to 




Figure 4.3. Linear Motion Cueing Filter Solution Procedure. 


4.2.2. Lateral Mode 

For the lateral (roll/sway) mode, the algorithm development is analogous to the 
longitudinal mode. In Eq. (4.1), the inputs 0 and a x are replaced by <p and a y 


respectively. The sensed rotational motion 9 in Eq. (4.2) is replaced by <j ) . The specific 
force f x and sensed specific force f x become / and / respectively, with Eq. (4.6) 
now computed as 

fy = a- g<p - R Sz 0. (4.26) 

These changes thus result in the differential equation given in Eq. (4.9) becoming 
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(4.27) 


Z.+K + bjA + W, 


GqtoKqtO 


R sz ( A) - B o - B M\ ~ (s + R sAA )«i - s4> \ U A + u 2 + A (] i 


which when rewritten as Eq. (4.10) , will produce the state equation representation for the 
otolith model similar to Eq. (4.11), with the same system matrices except 
®oto = \Goto K'qto r s z 0] . The state equation representation for the vestibular model of 
the form of Eq. (4.13) ultimately results. For this mode the additional platform states 

given in Eq. (4.14) are now x d = [JIM * 3 jja y dt 2 !«/// . 


The remaining development is identical in form to Eqs. (4.15) to (4.25), resulting 
in a matrix of seventh-order transfer functions W(s) for the lateral mode. The on-line 
implementation of this mode is identical to Figure 4.2. 

4.2.3. Vertical Mode 

For the vertical, or heave mode, the single degree-of-freedom input u = a z , with 
the specific force /, = a z -g. The otolith model given in Eq. (4.5) then becomes 


(s + 4>) 


fz ( 5 ) Goto R OTO l , D \t . n \ 

(s+B o p + R l ) 


fzi s )’ 


(4.28) 


and can then be defined in state space notation as 
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where x OTO are the otolith states for this mode, and 
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Since this mode consists of a single translational degree-of-freedom, the formulation does 


not include the semicircular canals model, therefore A v = A OT o. B v = B OTO , and C v 
Cqto- This results in a sensation model of the same form as Eq. (4.13): 


x e =A v x e +B v u s -B v u A 
e = C v x e . 


(4.30) 


Similar to the longitudinal mode, additional motion platform states are included in the 


state equations: 


=A d x d +B d u s , 


(4.31) 


where x d represents the additional motion platform states: 

x d = [JIM 3 jja z dt 2 \ a -dt 


and A d and B d now become 
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The aircraft input u A now consists of a single channel of filtered white noise with break 
frequency y and can be expressed as 


K = - r*n + 7™ 


(4.32) 


The state equations given in Eqs. (4.30), (4.31), and (4.32) can then be combined 
to form the desired system equation of the same form as Eq. (4.16), where y is the desired 

output, and x = [x e x d x n ] 1 represents the combined states. The remaining 
development is identical in form to Eqs. (4.16) to (4.25), resulting in a single fourth-order 


63 



transfer function W 22 for the vertical mode. The block diagram of the on-line 


implementation of this mode is shown in Figure 4.4. 


Pa 



Figure 4.4. Optimal Algorithm Implementation for Vertical Mode. 


4.2.4. Yaw Mode 


For the yaw mode, the single degree-of-freedom input is u = !// . The state space 


representation is the same as Eq. (4.4) with output \j/ replacing 9 : 


v — A v -L R II 

see sec sec sec 


^ = C scc x N3 + D scc u. 


(4.33) 


Since this mode consists of a single rotational input, the formulation does include the 
otolith model, therefore A v = A scc , B v = B scc , C v = C SC c> and D v = D scc . This results in a 
sensation model of the same form as Eq. (4.13): 


x e = A v x e + B v u s - B v u a 

C CyX e DyUg DyU^. 


(4.34) 


Similar to the longitudinal mode, additional motion platform states are included in the 
state equations 


=A d x d + B d u s , 


where x d represents the additional motion platform states x d = j y/dt y/ 


(4.35) 


, and A d and 


B d now become 
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The aircraft input u A now consists of a single channel of filtered white noise with break 
frequency y and can be expressed as 


K = - 

U A = X „ 


(4.36) 


The state equations given in Eqs. (4.34), (4.35), and (4.36) can then be combined 
to form the desired system equation of the same form as Eq. (4.16), where y is the desired 

output, and x = [x e x d x n ] 1 represents the combined states. The remaining 


development is identical in form to Eqs. (4.16) to (4.25), resulting in a single fourth-order 
transfer function W n for the yaw mode. The block diagram for the on-line 
implementation for this mode is shown in Figure 4.5. 



Figure 4.5. Optimal Algorithm Implementation for Yaw Mode. 


4.3. Pilot Tuning of the Algorithm 

A set of motion cueing filters for the longitudinal, lateral, vertical, and yaw modes 
was developed using the solution procedure given in Figure 4.3. The new semicircular 
canals model given in Eq. (3.8) was implemented, along with the Young-Meiry otolith 
model from Eq. (3.11) that was previously implemented [13], [9]. Filtered white noise 
break frequencies were initially set at 1 rad/sec for each degree-of-freedom. In some 
instances, the MATLAB™ error message “cannot order eigenvalues; spectrum too near 
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imaginary axis” resulted from the motion cueing filter development. This error occurs 
when the Hamiltonian matrix for the Riccati equation has eigenvalues on or very near the 
imaginary axis. To resolve this, a small change to any of the weights was needed to 
obtain a successful solution to the Riccati equation as suggested by Brogan [41]. 
Nonlinear scaling coefficients for each degree-of-freedom were based upon those chosen 
by Wu [13]. 

In order to determine the nonlinear gain coefficients for each degree-of-freedom 
that resulted in the most desired pilot performance, a trained simulator pilot executed a 
series of pilot controlled maneuvers with the optimal algorithm on the NASA Langley 
Visual Motion Simulator (VMS) described in Section 2.1. A series of maneuvers were 
first executed with the coefficients determined prior to testing. Coefficients for each 
degree-of-freedom were then adjusted until the simulator pilot subjectively felt the 
desired perception and performance were reached, while ensuring that the simulator 
motion platform limits were not exceeded. The following maneuvers were executed for 
the algorithm: 

Straight Approach and Landing (with varying wind from head to tail) 
Offset Approach and Landing (with and without turbulence) 

Pitch, Roll, and Yaw Doublets 
Throttle Increase and Decrease 
Coordinated Turn 

Ground Maneuvers (taxiing, effect of aircraft brakes) 

Takeoff from Full Stop. 

The optimal algorithm resulted in motion cues with which the simulator pilot 
commented he had more control and confidence as compared to the NASA adaptive 
algorithm. For both pitch and roll doublets, a fast response was observed when changing 
directions. On takeoffs, the optimal algorithm was found to be easier to pitch up to the 
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desired attitude and control the aircraft. A noticeably large side force was observed with 
the coordinated turn maneuver. By reducing the gains for the roll degree-of-freedom, this 
side force was reduced to a minimal sensation. The pitch gains were decreased to reduce 
the likelihood of entering the braking region or exceeding the actuator limits. Reducing 
the gains for both roll and pitch degrees-of-freedom still yielded desirable motion cues. 

The severe turbulence effects that were included with the offset approach and 
landing maneuver were hardly noticeable. An increase of the vertical gain coefficients 
resulted in increased cues, but still less than satisfactory. This increase in the vertical 
gains (coupled with an increase of the surge gains) resulted in forward surge cues that are 
more coordinated with the pitch cues, and a larger aft surge cue (initially, the aft cue was 
noticeably smaller than the forward cue). 

The effect of the otolith model upon the vertical filter characteristics was 
investigated. Figure 4.6 compares the frequency response of the original heave filter 
using the Young-Meiry model with the response using the proposed otolith model given 
in Eq. (3.29). Note that the original filter results in a gain decrease starting at about 5 
rad/sec, while the proposed model filter shows the gain unchanged for the same high 
frequencies. For the original filter, the filtered white noise break frequency y was 
increased to 471 rad/sec (2 Hz) to remove a right-half plane zero that resulted in a large 
false cue at both the onset and end of the pulse. For the revised filter with the proposed 
otolith model, this break frequency was reduced to 1 rad/sec, resulting in the specific 
force cue shown in Figure 4.7. The proposed model filter results in a faster onset cue that 
approaches the aircraft step input, and a faster washout that reduces the maximum 
simulator displacement while sustaining the peak onset magnitude. 
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Vertical Filter Frequency Response 
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Figure 4.6. Vertical Filter Frequency Responses with Young-Meiry and Proposed 
Otolith Models. 
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Figure 4.7. Specific Force Responses with Young-Meiry and Proposed Vertical 
Filters. 


Simulator Z-Displacement 
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A filter matrix W(s) is then generated for the pitch/surge mode using the proposed 
otolith model. In order to eliminate a singular response, both the lead time constant t l 
and the short time constant t 2 are removed from the semicircular canals model 
implemented in the motion cueing algorithm: 



(1 + v)(l+ v) 


u , . 


(4.37) 


Eq. (4.37) can be rewritten as 


9 = 


G scc s 


s~ + T,s + T n 


Ui , 


(4.38) 


where T 0 and T t become 


T - 1 

1 0 — 9 — 

Wi 1 


and can be defined in state space notation as 


*scc ^scc x scc ®scc u 


Q ^"SCC X SCC ®SCC U ’ 


(4.39) 


where in observer canonical form, 

’ C scc - [1 0]> and D scc — [G scc 0]. 

There is one less state as compared to Eq. (4.4), which in turn results in a matrix 
of sixth-order filters for W(s). The translational break frequency y 2 was increased from 1 
to 4k rad/sec for the original filters, eliminating a false specific force cue at onset. For 
the proposed model filters, y 2 is reduced to n rad/sec to produce a faster onset cue. The 
semicircular canals threshold G SC c was reduced to 2.0 deg/sec, the value obtained by 
Benson, et al. [37], to reduce the magnitude of the tilt coordination rate. Figure 4.8 
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compares the responses to a ramp to step input for both filters. Using the proposed 
otolith model in the filter development results in a faster onset ramp for the specific force 
response, with a faster onset and reduced magnitude for the tilt coordination rate. For 
both the original and revised filters, the weight component Q(2,2) needed to be increased 
from 1 to 10 to produce the magnitude of the specific force cues shown in Figure 4.8. 
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Figure 4.8. Responses to a Surge Ramp to Step Input with Original and Revised 
Longitudinal Filters. 
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Figure 4.9 compares responses to a lateral (sway) half sine input for both the 
original and revised filters. Since the algorithm development is similar to the 
longitudinal mode, the effect of the change of otolith models upon the motion cues is 
expected to be the same as the longitudinal mode. Note that the revised filter has a faster 
onset ramp with a larger specific force cue. 
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Figure 4.9. Responses to a Sway Half Sine Input with Original and Revised Lateral 
Filters. 


A summary of the optimal algorithm parameters for each of the four modes 
(longitudinal, lateral, vertical, and yaw) is given in Appendix B, Table B.l. The filter 
characteristics for both the original and revised linear motion cueing filters W(s) are also 
given in Appendix B, in Tables B.2 and B.3 respectively. 
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5. Integrated Human Perception Model 

The purpose of this study is to develop a model of human motion perception that 
can be readily implemented into the nonlinear motion cueing algorithm discussed in the 
next chapter. This integrated human perception model includes both vestibular and 
visual motion stimuli and incorporates the interaction between the vestibular and visual 
stimuli. This study is based on the literature presented by several researchers who 
investigated both the characteristics of visually induced motion sensation and the visual- 
vestibular interaction. 

5.1. Visually Induced Self-motion 

The visually induced effect on self-motion perception is commonly known as 

vection. Circularvection refers to visually induced rotational motion, in particular yaw, 

but also visually induced pitch and roll. Linearvection refers to visually induced 

translational motion. One common experience of linearvection is the illusion of moving 

backwards when seated in a stationary train car as the adjacent train in the station begins 

to slowly move forward. The self-motion response to a full visual field surround rotating 

about a vertical axis has been described by Young [42]: 

The response to a full field surround which suddenly begins to move at 
constant velocity is rather startling, although quite repeatable. At first, the 
veridical motion is sensed - the surround appears to be moving and the 
subject feels himself stationary. After a period of typically two to five 
seconds, the visual field appears to slow down, often to a stop, and the 
subject perceives himself as rotating in the opposite direction. The 
sensation of rotation builds to a maximum over a period of three to ten 
seconds, rising approximately as an exponential. 

Young [42] then noted that in order to achieve a complete “saturation” of this 
effect, in which the visual field is perceived to be entirely stationary, it is useful to have a 
wide, compelling field of view in the periphery, moving uniformly at speeds less than 60 
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degrees per second. Young then noted that if the visual surround is allowed to accelerate 
smoothly to its final velocity, at accelerations comparable to the acceleration thresholds 
of the semicircular canals, then the self-motion is more likely to be perceived as a 
smooth, continuous development of circularvection. 

Visually induced self-motion has been explored for rotations about both the earth- 
horizontal and earth-vertical axes, and along all three linear axes. The general 
characteristics of visually induced self-motion in the absence of confirming vestibular 
stimuli have been reported by Young [43] and supported by other researchers. Young 
noted two distinct classes of visual cues for flight simulation: the foveal cues, the high 
acuity, high information-dense central field cues that must be “read” to be interpreted, 
and the peripheral cues, the wide-field, lower acuity, rapidly moving cues that generate 
non-cognitive motion perception. These cues correspond respectively to the high static 
acuity, cone-filled fovea, and the high dynamic sensitivity, rod-filled periphery of the 
retina. 

Brandt, et al. [44] demonstrated that the peripheral visual field was of primary 
importance in stimulating self-motion over the central visual field. They observed that 
when the central visual field is masked up to 120 degrees in diameter, circularvection 
diminished very little. Conversely, if peripheral vision was precluded, stimulation of the 
central field of up to 60 degrees in diameter produced an almost exclusive exocentric 
perception of the moving visual surround and a very weak self-motion perception. They 
also found that when equal stimulus areas are presented either foveally or peripherally, 
that stimulation of the peripheral visual field is more favorable to stimulate self-motion 
perception. 
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Brandt, et al. [45] demonstrated that for circularvection, background stimulation 
is dominant over foreground stimulation. They showed that movement in the background 
induced the apparent self-motion, while if the foreground moved the stationary 
background inhibited circularvection. They showed that if a stationary pattern is attached 
to the background, the circularvection latencies are increased significantly, thus 
indicating an inhibitory effect of circularvection due to the presence of stationary 
contrasts in the background. Howard and Howard [46] later demonstrated that the 
presence of stationary objects in the foreground significantly increases circularvection 
and reduces the latency of onset to circularvection. 

The spatial frequency of the scene determines its effectiveness in generating self- 
motion. Held, et al. [47] demonstrated this by quantifying the visual border placement 
and velocity necessary to achieve a visually induced roll. In their experiment the 
observer viewed a circular disc through a monocular color-corrected lens, effectively 
producing an extended visual field. The disc was covered with a random pattern of spots 
with areas that were masked off, producing a set of ring-shaped displays. They found 
that, in general, the magnitude of visually induced tilt increased with field size, with the 
use of a large number of rings subtending small solid angles in the peripheral areas being 
more effective than the same for central areas. Young [43] commented that the 
peripheral field display should also have a sufficient number of borders such as stars, 
clouds, or ground features to induce the perceived self-motion. 

The visual field velocity determines the magnitude of the self-motion up to a 
saturation velocity that most likely corresponds to the blurring of the visual field 
associated with increased dynamic acuity [43]. Saturation of vection occurs when the 
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field appears stationary in space and all motion is perceived self-motion or egocentric. 
Brandt, et al. [45] observed that the velocity of apparent yaw self-motion matches the 
stimulus speeds up to 90 to 120 degrees per second. Young and Oman [48] observed for 
a pitch stimulus, the visually induced pitch increases with field speed up to a stimulus of 
about 40 to 60 degrees per second. They also observed a limit of about 40 to 60 degrees 
per second for visually induced roll that was confirmed by Held, et al. [47]. Young [43] 
noted that saturation occurs for linearvection up to 1 meter per second, which was the 
maximum stimulus velocity tested Berthoz, et al., [49]. 

Young [43] found that the approximate frequency response for both 
circularvection and linearvection is flat from static inputs up to a frequency of 0.1 Hz, 
beyond which it decreases at least as rapidly as a first-order filter. Berthoz, et al. [49] 
confirmed these results for forward linearvection, with similar results obtained for 
visually induced pitch [50] and for yaw circularvection [51]. 

5.2. Latency of Onset to Yection 

Young [43] noted that the onset delay, or the latency, of visually induced motion 
is highly variable among individuals. Repeated exposures will reduce this latency, as 
does the development of the appropriate mental set, thus allowing for the development of 
vection. The latency of onset to either circularvection or linearvection has an impact on 
the perception of motion in flight simulation. Several experimenters have quantified this 
phenomenon. 

Brandt, et al. [44] conducted experiments using a rotating chair located in the 
center of a closed cylindrical drum 1.5 meters in diameter, whose inner walls were 
painted with alternating vertical black and white stripes subtending 7 degrees of visual 
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angle. To stabilize the direction of the chair, the subjects were asked to fixate on a one- 
degree luminous spot mounted on the chair and presented in a position straight ahead of 
the subject. The yaw stimuli consisted of drum rotations moving at either constant 
velocity (10 to 180 deg/sec) or at constant acceleration (1 deg/sec 2 ). The latency of 
onset was measured using a stopwatch from the sudden onset of optokinetic stimulation 
(lights on) to the beginning of circularvection. Following the stimulus onset, 
circularvection began after a latency of about 3 to 4 seconds. The latencies were found to 
be independent of the stimulus velocities tested. 

Young and Oman [48] carried out experiments in one of the differential 
maneuvering simulators (DMS) at the NASA Langley Research Center. Each simulator 
consisted of a jet cockpit mounted on a fixed-base platform inside of a forty-foot 
diameter projection sphere. Visual scenes were projected on the interior of the sphere 
wall by a computer controlled projection system that consisted of two servo-driven 
transparent plastic hemispheres on which the scenes to be projected were painted. A high 
intensity point light source mounted near the center of each hemisphere projected the 
scene onto the interior walls of the simulator sphere. The hemispheres projected a pattern 
of randomly spaced and oriented black rectangles of 2 to 3 degrees in subtended angle 
against a white background, with a black-white ratio of approximately 25%. A series of 
constant velocity yaw stimuli were presented randomly left and right at speeds of 5 to 
120 deg/sec. The latency of onset was recorded using a stopwatch. A rapid decrease in 
time to onset of circularvection with increasing pattern speed (from 11 seconds at 5 
deg/sec to 6 to 2 seconds from 10 to 120 deg/sec) was observed. 
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Howard and Howard [46] performed tests using an apparatus consisting of a 
vertical cylinder of translucent plastic, with an inner surface covered with adhesive black 
vinyl and perforated with round holes randomly distributed over the surface. The 
cylinder was illuminated from outside by diffuse tungsten light. The visual field of the 
subject inside the cylinder was filled with a random array of white spots subtending 
approximately 0.5 and 1 degree respectively. The cylinder rotated from left to right, from 
the subject’s point of view, about its mid- vertical axis. Six stimulus field conditions were 
tested at angular velocities of 5, 10, and 25 deg/sec. The field conditions consisted of a 
full field condition where subjects saw only a moving display without stationary objects, 
a set of conditions with two vertical rods placed symmetrically at central, intermediate, 
and peripheral locations in front of the moving display (each presented separately and all 
presented together), and a condition with a frame placed in front of the moving display 
similar to the “window bars” that frame the video monitors in a simulator cockpit. Each 
field condition at each velocity was tested once with the subject looking straight ahead 
with relaxed gaze, and once with the gaze fixated on a stationary white spot projected 
from a laser and positioned at eye level straight ahead of the subject. The latency of 
vection was measured by having the subject press a switch at the first sign of vection. 

The results obtained by Howard and Howard [46] showed that latency is longer 
when there are no stationary objects in view. They note that this effect is most evident at 
the lowest stimulus velocity, where subjects were usually unaware that the display was 
moving. They reasoned that at this velocity the eyes reflexively pursue a moving display 
without the presence of stationary objects. Fixation upon a small laser point was also 
sufficient to increase the vection magnitude significantly. At 5 deg/sec, the presence of 
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the visual frame reduced latency for the full field condition from 48.4 seconds to 8 
seconds with no fixation. With fixation, the presence of the visual frame reduced latency 
from 14.3 seconds to 5.4 seconds, with the latency relatively unchanged for this condition 
(5.2 sec) for 10 and 25 deg/sec. The condition of central vertical rods without fixation 
also yielded results within the ranges reported by Young and Oman [48]; a latency of 9.4 
seconds was found for a stimulus of 5 deg/sec that decreased to 5.6 seconds for 25 
deg/sec. 

Berthoz, et al. [49] tested the latency to onset of forward linearvection. The 
experimental apparatus they used allowed the projection of a moving 35-mm film loop of 
randomly distributed images on a screen that was fixed on a mobile cart. The screen 
image was projected via two mirrors to produce two peripheral images parallel to the 
sagittal plane of the head. The subject, whose head position was fixed by a chin rest, 
could view the moving images through a black box with side windows that limited the 
visual field between 20 and 70 degrees away from the sagittal plane on each side. The 
sensation of self-motion experienced by the subject was measured by the method of 
magnitude estimation by adjusting a lever fixed to the cart that could rotate forward or 
backward starting from a zero position. Both the lever rotation and the image velocity 
were recorded with a potentiometer. 

In the experiment, latencies of about 1 to 1.5 seconds were observed for velocities 
measured between 0.2 and 1 m/sec. This significant difference in latency between 
linearvection and circularvection may be related to the differences in response dynamics 
associated with the otolith and semicircular canals respectively. 
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5.3. Visual-Vestibular Interaction Models 


Zacharias [52] reported that both psychophysical and neurophysiological studies 
support the theory that visual and vestibular cues are jointly processed to provide for a 
perceived sense of self-motion. The vestibular nucleus complex was identified as a 
possible interaction for the convergence of sensory inputs. Zacharias [52] then noted that 
experiments reported by Dichgans, et al. on the single unit recordings from the vestibular 
nuclei of goldfish indicated that the majority of cells respond to both vestibular and 
moving visual field inputs. When both visual and vestibular stimuli were presented in 
opposing directions consistent with rotation in the presence of a physically stationary 
visual surround, the afferent firing rate was characterized by the faster response and 
greater sensitivity of vestibular stimulation combined with the non-adapting behavior of 
visual stimulation. The result was a signal that accurately indicated the perceived angular 
velocity. 

A study by Young, et al. was also reported by Zacharias [52] in which subjective 
velocity and acceleration were measured in response to combined yaw-axis rotational 
cues. The study showed that the subjective velocity response was biased in the direction 
of the induced circularvection, but not to the extent of a simple summation of 
circularvection and expected vestibular response. These studies indicated that a simple 
linear summation of the visual and vestibular cues failed to predict the response when 
both cues are simultaneously presented. 

Visual motion cues dominate the perception of velocity and orientation in the 
steady state and at frequencies below 0.1 Hz [43]. At higher frequencies, the vestibular 
cues will tend to dominate. Confirming vestibular cues, in the direction opposite to the 
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visual field, can produce a rapid onset of the visual self-motion that is sustained by vision 
after the vestibular cues have been washed out. When visual and vestibular motion cues 
are in conflict, either due to the direction of motion or to a difference in magnitude, the 
vestibular cues will initially dominate. Young [43] suggested that when both inputs are 
presented to a subject simultaneously, he or she would combine or “mix” the two cues in 
a nonlinear manner, favoring the visual input for confirming cues and the vestibular input 
for conflicting cues. 

Young [53] first proposed that visual and vestibular cues are independently 
processed to produce two estimates of motion that are compared with one another to 
provide some means of cue conflict. For low conflict, such as when the cues confirm one 
another, the perceived motion is calculated from the weighted sum of the two estimates. 
This weighting is dependent on the sensory cue characteristics in the given situation and 
would be chosen to minimize the error in the combined cue estimate. For high conflict, 
that is when the cues fail to confirm one another, the weighting is then shifted based upon 
the reliability of each cue. 

Zacharias [52] developed a cue conflict model for yaw perception that was based 
on the switching concept first proposed by Young [53]. This model is illustrated in 
Figure 5.1. For low conflict, that is when the visual and vestibular cues are confirming, 
the perceived motion is calculated from a weighted sum of the two estimates. For high 
conflict, that is when the cues fail to confirm one another, the weighting on the visual 
input is reduced and that on the vestibular input is increased until the conflict is reduced. 
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Figure 5.1. Visual-Vestibular Conflict Estimator Model. (Zacharias, [52]) 

The visual cue is passed through an internal model of the vestibular dynamics (a 
reduced model with only the long time constant was used) to produce an “expected” 
vestibular signal that is then subtracted from the actual vestibular signal. To allow for 
long-term resolution of steady-state conflict, the absolute value of this error is passed 
through an adaptation operator to generate the conflict signal. The adaptation operator 
determines how long a steady state conflict is allowed to continue by washing out the 
conflict signal, ultimately allowing for an averaging of the two cues. Zacharias [52] 
chose a value of 10 seconds for the adaptation time constant based on typical trainer 
acceleration latencies observed for a conflicting visual field. 

The weighting of each cue is governed by a gain K that is derived directly from 
the two cues, and varies between zero and one. The gain K is computed from symmetric 
weighting functions that are applied to the vestibular and the washed out conflict signal 
Zacharias [52] noted that “for simplicity” a cosine bell operator was chosen for the 
weighting function. A large conflict signal will drive the gain K to zero, gating out the 
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visual path, while a small conflict signal will drive K to a peak weighting value between 
0.5 and 1, depending on the amplitude of the vestibular signal. In a low conflict situation 
the cues are either averaged or the visual cue is passed straight through and the vestibular 
cue is gated, depending on the magnitude of the vestibular cue. Zacharias [52] then 
noted that the value of the conflict threshold e was determined by appealing to past work 
in defining threshold behavior. He presumed that the same type of behavior associated 
with the vestibular motion thresholds characterized cue conflict detection. 

A second approach to modeling the visual-vestibular interaction was developed 
by Borah, et al. [54]. This approach involves the implementation of an optimal estimator 
as a “central processor” representation for the central nervous system processing of 
sensory inputs. Individual sensory dynamic models represent the visual and vestibular 
systems, and this concept can be extended to include proprioceptive and tactile models. 
The sensors respond to input stimuli and send signals to a central processor represented 
by a steady-state Kalman filter, which combines the sensory information to generate an 
estimate of the perceived motion. In this model, a modified version of the cue conflict 
estimator proposed by Zacharias [52] was also implemented. 

Van der Steen [51] proposed a self-motion perception model in which vestibular 
and visual stimuli are combined to describe perceived self-motion. This model is shown 
in Figure 5.2. The model can describe perceived self-motion induced by either vestibular 
or visual stimuli alone, or a combination of both. However, unlike the model proposed 
by Zacharias [52], cue conflict estimation is not considered in this model. 
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Figure 5.2. Model for Self-Motion Perception with Visual Attractor. (Van der 
Steen, [51]) 

Van der Steen [51] introduced the concept of a neural filter in the model. The 
neural filter transfers the afferent response of either the visual or vestibular sensor to a 
perceptual physical variable. In other words, the neural filter models the process of 
assigning a perceptual meaning to a sensory output signal. Van der Steen then noted that, 
in general, neural filters are not explicitly described in the literature but are either 
represented as a constant that relates the afferent response to the perceived response, or 
are implicitly imbedded in the sensory dynamics. For example, the vestibular transfer 
function H VEST cascaded with the vestibular neural filter NF VEST represents the perceived 
self-motion estimate from vestibular stimuli. 

Van der Steen [51] then noted that psychophysical experiments concerning 
vection showed that, depending on self-motion that the visual scene suggests, the visual 
estimate of self-motion “attracts” the vestibular estimate of self-motion. He suggested 
that the self-motion perception model needs a component that handles this “attraction” 
towards the visual system’s estimate of self-motion. This component was referred to as 
the “visual attractor”. 
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The visual attractor uses the visual and vestibular system’s estimates of perceived 
self-motion. The difference between these cues is then filtered as shown in Figure 5.2, 
forming an optokinetic influence that is an estimate of perceived self-motion from visual 
stimuli. Van der Steen [51] noted that the filter H OK represents the gradual build-up of 
perceived self-velocity when exposed to a constant visual scene velocity as observed in 
psychophysical experiments, and can be represented by a first-order low-pass filter of the 
form 



The perceived self-motion yielded by the model is then the sum of the optokinetic 
influence and the vestibular system’s estimates of perceived self-motion. 

Van der Steen [51] determined values for the optokinetic time constant r OK from 
experiments conducted using an optokinetic drum. The apparatus consisted of a rotating 
chair surrounded by a closed cylindrical drum, the inside of which contained alternating 
black and white stripes. Each test subject was asked to fixate on a stripe edge near the 
middle of the drum to indicate left or right drum motion. Two experiments were 
performed using this apparatus. In experiment 1, six visual acceleration amplitudes were 
tested with the chair stationary, and fourteen inertial acceleration amplitudes were 
provided with three constant magnitude drum accelerations. In experiment 2, the drum 
was accelerated for one second to a constant velocity of either 10 or 20 deg/sec, with the 
chair remaining stationary. After 17 seconds, the drum decelerated for three seconds. 
Four acceleration amplitudes were tested. In both experiments, the subject indicated 
perceived drum motion by pushing a button, with the elapsed time recorded 
electronically. 
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Van der Steen [51] then performed a series of simulations of the self-motion 
perception model with the data from both experiments. It was found that the model could 
be more accurately described with different optokinetic time constants for each 
experiment; the value for t ok was chosen smaller for experiment 1 (2 seconds) and larger 
for experiment 2 (10 seconds). A smaller time constant for smaller drum accelerations 
also resulted. Van der Steen concluded that the model dynamics using fixed parameters 
might not completely describe results for each combination of stimuli. 

5.4. Visual Sensory Dynamics Models 

Each visual-vestibular interaction model examined incorporates a model of the 
visual receptor dynamics that in turn produces a perceptual estimate of the visual scene 
motion. Zacharias [52] did not model visual sensory dynamics due to the lack of 
experimental data for single channel visual response, and assumed that the visual system 
has a relatively wide-band response. The negative sign in Figure 5.1 reflects the fact that 
the visual field motion is opposite in direction to the perceived self-motion, i.e., a visual 
field moving to the left induces self-motion of the subject to the right. 

Borah, et al. [54] modeled the dynamics of the visual sensor as unity, noting that 
the eye detects the visual field motion almost immediately after a short neural 
transmission delay. Van der Steen [51] modeled the perceptual dynamics as a cascade of 
the visual receptor transfer function and neural filter with a unity gain and a delay y. 

NF VIS H VIS =-e^. (5.2) 

Hosman and Van der Vaart [55] noted that r d is due to the delay of the visual 
receptors along with the delays due to both neural transmissions from the retina to the 
visual cortex and information processing during motion perception. From experiments in 
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roll rate perception with visual displays only, Hosman and Van der Vaart found values 
for T d to be about 90 msec for peripheral visual field stimulation and about 150 msec for 
central visual field stimulation. 

5.5. Proposed Rotational Model 

A revised visual-vestibular interaction model will now be constructed for 
rotational motion. This model can be used to estimate perceived motion for yaw, roll, 
and pitch stimuli. As suggested by Borah, et al. [54], the visual motion cues considered 
will be limited to peripheral visual scenes provided by a flight simulator with a wide 
visual scene field. These peripheral cues would be equivalent to the passage of stars or 
clouds in a wide field simulation. The cues do not include any elements in the structure 
of the scene such as landmarks, orientation cues, or a visual horizon. 

A visual-vestibular interaction model for rotational motion is proposed and is 
shown in Figure 5.3. The proposed semicircular canals model given in Eq. (3.9) is used. 
The vestibular model combines the afferent dynamics model with the neural filter gain 
proposed by Van der Steen [51], resulting in a model with a perceived response to 
vestibular stimuli. Since the visual motion cues are assumed to be peripheral, the visual 
delay T d = 90 msec obtained by Hosman and Van der Vaart [55] is chosen. The 
optokinetic influence proposed by Van der Steen [51] is also implemented. 
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Figure 5.3. Proposed Visual-Vestibular Interaction Model for Rotational Motion. 

As first proposed by Zacharias [52], the model produces a washed out conflict 
signal co err . The weighting of the optokinetic gain K ok is then computed by a modified 
cosine bell function suggested by Borah, et al. [54] as shown in Figure 5.4. The gain K ok 
varies between zero and one. A conflict signal greater than the threshold value £ {co err > 
£), will drive the optokinetic gain to zero, whereas a signal below the threshold value 
{(Oerr < £) will drive the gain to a value between zero and one, approaching one as co err 
approaches zero. For co err < 0, the gain remains at one. As previously suggested by 
Borah, et al. [54], the vestibular path gain remains fixed at unity. 

The conflict threshold £ is chosen to equal the vestibular indifference motion 
threshold [52]. The angular velocity thresholds obtained by Benson, et al. [37], 1.6 
deg/sec for yaw stimuli, and 2.0 deg/sec for roll and pitch stimuli, will be used in the 
model. 
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Figure 5.4. Modified Cosine Bell Operator for Optokinetic Gain. 

In order to examine model responses to various stimuli, a MATLAB/SIMULINK 
representation of the proposed rotational model shown in Figure 5.3 was developed. 
Time constants t ok = 2 seconds and r r = 5 seconds were chosen to produce the latencies 
noted in the literature. In the model, the latency is defined as the amount of time to 
perceive motion above a visual indifference threshold of 3 deg/sec. Model responses to 
yaw inputs with either visual cues alone or confirming visual and vestibular cues were 
examined. 

Figure 5.5 shows the responses to a visual field step input of 10 deg/sec. Since 
there is no vestibular input, the rectified error is the magnitude of the visual input filtered 
through the internal model of the semicircular canals. The adaptation operator then 
generates the washout error co err . Due to the large value of co err , the cosine bell function 
will produce a gain of zero for nearly five seconds. This results in a corresponding 
latency in the perceived angular velocity response. 
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Figure 5.5. Rotational Perception Model Responses to Visual Step Input of 10 
deg/sec. 

Once co en decreases below the conflict threshold e ; the optokinetic gain will vary 
between zero and one, resulting in the onset of perceived motion or circularvection. This 
gain will rapidly rise to a value of one once co en reaches zero. As (O err becomes negative, 
the gain K ok remains at one. If a cosine bell operator were applied to this negative 
response, the gain would decrease back to zero, resulting in a large sag in the perceived 
response. The perceived motion reaches its maximum value with a rise time of about ten 
seconds, as governed by the time constant z OK . 

Various magnitudes of angular velocity inputs were examined in order to compare 
latency responses with those obtained from psychophysical experiments in the literature. 
Figure 5.6 compares the model responses to visual step inputs of 5, 10, and 25 deg/sec. 
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Rotational Model Response to Visual Step Inputs 
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Figure 5.6. Rotational Model Responses to Visual Step Inputs of 5, 10, and 25 
deg/sec. 

Due to the conflict estimator, the responses in Figure 5.6 produce a “dead zone”, 
of which the duration increases with increasing step input. Assuming this phenomenon 
as latency alone contradicts the experimental results obtained by Howard and Howard 
[48] and Young and Oman [46] that showed the latency decreases with increasing stimuli 
magnitudes. However, assuming latency to onset occurs until an indifference threshold is 
reached reveals the latency decreasing with increasing step inputs. 

The latencies resulting for the proposed rotational model are shown in Table 5.1. 
The latencies obtained with a threshold of 3 deg/sec result in values that are near those 
obtained by Howard and Howard [46] (5.2 to 5.4 seconds) for the stationary visual frame 
condition with fixation. As seen in Figure 5.6, increasing the visual threshold to 4 
deg/sec for a 5 deg/sec step input would result in a latency of 7 seconds, which 
approaches the value of 8 seconds obtained by Howard and Howard [46] for the same 
condition with no fixation. 
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Table 5.1. Model Results for Latency to Onset of Circularvection. 



Latency (sec) 

(this (deg/sec) 

Model 

Howard and Howard 

5 

5.525 

5.4 

10 

5.0 

5.2 

25 

5.0 

5.2 


Figure 5.7 shows the model responses to a confirming visual and vestibular step 
inputs of 10 deg/sec. Due to the visual delay r d , a large value of co err results at the onset 
that is rapidly washed out in a fraction of a second, resulting in the optokinetic gain K ok 
increasing from zero to one during this instant and remaining at one for the duration of 
the response. Due to the rapid onset of the semicircular canals, the visual delay will have 
a negligible effect on the perceived response. After the vestibular cue decays, the 
optokinetic influence will gradually increase until the maximum response is achieved. 
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Figure 5.7. Rotational Perception Model Responses to Confirming Step Inputs of 10 
deg/sec. 
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5.6. Proposed Translational Model 

A visual-vestibular interaction model was developed to estimate motion in all 
three translational axes. The same assumptions applied to the rotational cues were 
considered. The proposed visual-vestibular motion model for translational motion is 
shown in Figure 5.8. The model structure is similar to the rotational model. The same 
values for visual delay T d and optokinetic time constant z OK as proposed for the rotational 
model will be used. In this model the washout error v err is used to estimate the 
optokinetic gain K ok . Vestibular and optokinetic responses are combined to produce 
perceived linear velocity. 
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Figure 5.8. Proposed Visual-Vestibular Interaction Model for Translational 
Motion. 


A MATLAB™/SIMULINK™ representation of the model shown in Figure 5.8 
was developed. Figure 5.9 shows responses to a visual field step input of 1 m/sec. An 
adaptation time constant T c = 0.2 seconds was chosen to generate latencies close to those 
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obtained by Berthoz, et al. [49]. The rectified error is the magnitude of the visual 
velocity response filtered through the internal model of the otoliths. As a result of this 
fast time constant, the washout error decays very quickly, resulting in a small latency of 
about 1.5 seconds. The perceived linear velocity then reaches its maximum value in 
about ten seconds, as governed by the time constant t ok - 
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Figure 5.9. Translational Perception Model Responses to Visual Field Step Input of 
1 m/sec. 

Various magnitudes of linear velocity inputs were examined in order to compare 
latency responses with those obtained from psychophysical experiments in the literature. 
The latencies for T c = 0.2 seconds and s= 0.2 m/sec result in values that fall within the 
range (1 to 2 m/sec) reported by Berthoz, et al. [49] for velocity inputs from 0.4 m/sec to 
1 m/sec. As the input is reduced towards the 0.2 m/sec threshold, the latency increases 
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beyond the reported experimental range. Table 5.2 lists the latencies the model generates 


for inputs from 0.2 to 1 m/sec. 


Table 5.2. Model Results for Latency to Onset of Linearvection. 


v vis (m/sec) 

Latency (sec) 

0.4 

1.75 

0.6 

1.25 

0.8 

1.05 

1.0 

0.975 


Figure 5.10 shows the model responses to confirming visual and vestibular step 
pulse inputs of 1 m/s 2 magnitude and 1- second duration, which produces a ramp to step 
velocity input of 1 m/s. Due to the visual delay T d , a sub-threshold value of co err is 
generated at the onset that will result in the optokinetic gain K ok being slightly less than 
one for about one second. 
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Figure 5.10. Translational Perception Model Responses to 100% Confirming Pulse 
Inputs of 1 m/sec Magnitude and 1 second Duration. 
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Due to the rapid onset of the otoliths, the visual delay will have a negligible effect 
on the perceived response. After the initial otolith onset decays to its steady state value, 
the optokinetic influence then gradually increases until the maximum response is 
achieved. 
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6. Nonlinear Motion Cueing Algorithm 

6.1. Problem Description 

The motion perception process depends on the interaction between visual and 
vestibular sensation. Based on this, the problem is to develop a set of motion cueing 
filters that minimize the pilot perceptual error. A nonlinear approach is desired to further 
maximize the available motion cues within the hardware limitations of the motion 
system. This algorithm must be efficient enough to update filter characteristics in real 
time. Cardullo and Kosut [56] suggested this approach, and Ish-Shalom [57] also 
proposed a similar algorithm structure. The structure of the proposed solution is 
illustrated in Figure 6.1. 


Aircraft 

States u. 




Figure 6.1. Proposed Solution for Nonlinear Motion Cueing Algorithm. 

The nonlinear algorithm incorporates models of the human vestibular sensation 
system, the new semicircular canals and otolith models, along with the new integrated 
visual- vestibular perception model. A nonlinear control law is implemented to generate a 
scalar coefficient or that is a function of the motion platform states. The matrix Riccati 
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equation is then solved in real time as a function of a. resulting in the feedback matrix 
needed to compute the desired motion cues. This results in nonlinear filter characteristics 
that sustain small motion cues for longer durations, and generates faster washout of large 
platform motions. 

6.2. Algorithm Development 
6.2.1. Longitudinal Mode 

The algorithm development for the longitudinal (pitch/surge) mode is given 
below. The input u is the same as given in Eq. (4.1) for the linear algorithm 
development: 


" e~ 


Mj 

_ a ._ 


u 2 _ 


( 6 . 1 ) 


where 0 is angular velocity, and a x is the translational acceleration, with each term 
respectively set equal to m, and u 2 . The reduced-order semicircular canals model given 
in Eq. (4.38) is used: 



G$cc s 

s~ + T,s + T n 


(6.2) 


where T 0 and 7) become 


T - 1 ±- 


l 0 — ’ -M 

T a T l 




where T a and r, are the same as those used in the optimal algorithm and given in the 

semicircular canals model of Eq. (3.8). G scc is the angular velocity threshold that scales 
the response to threshold units, using the threshold of 2.0 deg/sec obtained by Benson, et 
al. [37]. Eq. (6.2) is defined in state space notation as 
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(6.3) 


x scc — ^scc x scc ®scc u 

@ — ^"SCC X SCC ®SCC U ’ 

where in observer canonical form, 

’ ^scc — [l 0], and D scc — [G scc 0]. 


A - 

^scc 


T 1 
-T 0 0 


®scc _ 


T'scc F 6 
—Tr SC c? o 0 


In an attempt to produce the desired motion cues that most closely represent the 
perceptual behavior of the aircraft pilot, the confirming case of the integrated perceptual 
model (neglecting conflict estimation) will be incorporated into the perceptual channel. 
The visual delay was also neglected since it has only a small effect on the perceptual 
response. 

For a simulator pilot, the perceptual input u = u s and for the aircraft pilot is u = 
u A . Therefore, u = u s - u A can be considered as input to the cueing algorithm. For the 
perceptual error states x e = x s - x A , the input to the optokinetic influence of the integrated 
perception model given in Eq. (5.1) becomes 

yv / yv \ / yv \ yv yv 

9 L =[ 9 A -e s y[e A -e A ye A - 9 s 

~ i X \A + GsCC U \A ) — (''"IS + G S cc U \S ) ( 6 . 4 ) 

— ^SCC^IA ~ X l ~ Tl scC^lS • 

The output of the optokinetic influence given in Eq. (5.1) (for gain K ok = 1) becomes 


&OK =■ 


1 -0 


^ OK S 1 s + T 2 


(6.5) 


where T 2 = 1 jr OK . This can be defined as an additional state space term x 3 : 


x 3 — T 2 x 3 + T 2 6 ok 

— — T 2 x 3 + T 2 (G scc u m — x l — G scc u s , ) (6.6) 

— — T 2 x l — T 2 x 3 — G scc T 2 u 1 . 
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Eq. (6.3) now becomes 


x scc — ^scc x scc ®scc u 

@PE — Cscc x scc + D scc u, 


(6.7) 


with 




1 

0 ' 


-Gsccf 

0" 




A - 

^scc 

~T 0 

0 

0 

’ ®scc 

-GscJo 

0 

’ ^-SCC [l 

0 1 

]> ®scc 


-t 2 

0 

-T 

1 2J 


-g scc t 2 

0 





[Gsec 0], 


The sensed velocity v x is related to the stimulus specific force f x by the otolith 
model of Eq. (3.29): 

K(s) = G„ 0 K'oto t * b! /,(»)■ < 6 - 8 > 

s(s + B,) 

where the break frequency B x is neglected, and G 0T0 is the linear velocity threshold that 
scales the response to threshold units. As with the linear algorithm, the specific force is 

f x = a ,+ 8° + R s z &- (6-9) 


In terms of u x and u 2 , Eq. (6.9) is transformed into the Laplace domain, where 


f x {s) = u 2 (s) + \g-~ R Sz s Imj ( 5 ). 


(6.10) 


V ^ J 

Note that the form of the transfer function of Eq. (6.8) is similar to the form of the otolith 
model given in Eq. (4.5), with v replacing / as output, and the integrator replacing the 
short time constant term (5 + B x ) . Thus, Eq. (4.10) from the linear algorithm becomes 


E + B (f x = 

GqtO^OTO 


R Sz ( B 0 — Af)u x + gu x + gA (j | u x dt + u 2 + A 0 u 2 


(6.11) 
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and can be rewritten as v x + av-cu l + du t + e J u i dt + fu 2 + gu 2 , and is then defined in 


state space notation as 


X OTO — Aoto X OTO ®OTO U 

v x — C OTO x OTO + D oxo u, 


( 6 . 12 ) 


0 10 0 
0 -a 1 0 

A oto = 0 0 0 0 

0 0 0 0 

0 0 0 0 

C OTO = [1 0 0 10], 


d - ac 


t OTO T Sz 


f 

h - af 

0], 


The output of the optokinetic influence (for gain K ok equal to 1) becomes 


V = V = =— = — . 

^ OK S 1 s +T 2 


(6.13) 


where T 2 = \/t ok . This can be defined as an additional state space term x 9 : 


Xg T 2 Xg + T 2 ( X 4 X 7 + GqjqK OT qKUi ) 

— — T 2 x 4 — T 2 x 7 — T 2 Xg + T 2 G OTO K OTO r ; u x . 


(6.14) 


Eq. (6.12) now becomes 


X OTO AoTO X OTO ®OTO U 

v xpe ~ C OTO x OTO + D oto u, 


(6.15) 


0 1 0 0 0 0 

0 -a 1 0 0 0 

0 0 0 0 0 0 

0 0 0 0 1 0 

0 0 0 0 -a 0 

-71 0 0 -71 0 -r. 


C OTO = [1 0 0 10 1], 


d -ac 


T G r 

_ l^OTO'Sz 


[~ G oto K o 


f 

h-af 

0 


OTO r Sz 
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The representations in Eqs. (6.7) and (6.15) can be combined to form a single 
representation for the human perceptual model: 


x v = A v x v + B v u 

y PE _ CyXy+DyU, 


(6.16) 


where x v and y PE are, respectively, the combined states and perceived responses, and 
Ay, By, Cy, and D v represent the perceptual models as one set of state equations: 


A srr 0 

Ay = SCC 

v ft A 

V ^OTO 


, By = 


n ’V 0 C* 

°OTOj L V ^OTO 


, Dy = 


It is assumed that the same perceptual model can be applied to both the pilot in 
the aircraft and the pilot in the simulator as shown in Figure 6.1. We define the pilot 
perceptual error e, resulting in 


X e =AyX e +ByU S -ByU A 
e — CyX e + DyUg — DyU A . 


(6.17) 


Additional motion platform states and filtered white noise states defined in the 
linear algorithm development in Eqs. (4.14) and (4.15) are again used. The desired 
system equation given in Eq. (4.16) is then formed, with the cost function J given in Eq. 
(4.17). The system equation and cost function are then transformed to the standard 
optimal control form given in Eq. (4.18). The cost function is augmented with an 
additional term e 2al proposed by Anderson and Moore [58]: 


J' =e\[ e 2at (x T R'x + u ,I R 2 u')jt|, 


(6.18) 


where R' is positive definite, R 2 is positive semi-definite, and the scalar coefficient a 


represents a minimum degree of stability in the closed-loop system where a > 0 . 
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and the Riccati equation solution P(or) can be partitioned as 

P 11 (or) P 12 («) P 13 («) 

P(or)= P 21 (or) P 22 («) P 23 («) , (6-23) 

_ P 3l(«) p 23 («) p 33 («)_ 

where the partitions correspond to the partitions of the system matrix A. Reid and Nahon 

[9] noted that when computing K, only a subset of the elements of P is needed. 
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Substituting Eq. (6.23) and the expression for R , 2 from Eq. (4.18) into the expression for 


K given in Eq. (6.20) results in 

Kj {a) = R ' 1 [ByP u + B d r P 21 + DyQC v ] 

K 2 {a) = R , 1 [ByP 12 + B d r P 22 ] (6.24) 

K 3 (or) = R , 1 [ByP 13 + B d P 23 - DyQCy ], 

where, by symmetry, P 12 = P 21 . 

Noting that x n = u A , remove the states corresponding to the x n partition: 


X e 

i 


0 


0 

A, 


-B* 

0 


X e 

X d 

U. 


+ 


B v 

B, 


u s . 


(6.25) 


and substituting Eq. (6.22) into Eq. (6.25) results in 


X e 


Ay -ByK^flf) -ByK 2 («) 

X e 

4- 

'-B y (I + K,(a))' 

_ X d_ 


■B d Kj (a) A d - B d K, (a) 

_v 

i 

-B a K,(«) 


u 


A- 


(6.26) 


A nonlinear control law is chosen to make a dependent upon the system states: 

a= x dQ 2 x d> (6-27) 

where q 2 is a weighting matrix that is at least positive semi-definite. As the system 
states increase in magnitude, i.e. with large commanded platform displacements and 
velocities, then a increases, resulting in faster control action to quickly wash out the 
platform to its neutral state. For small commands there will be limited control action, 
resulting in motion cues being sustained for longer durations. The feedback matrix 
K ( a ) is then determined by solving the Riccati equation of Eq. (6.19) in real time as a 


function of a. 
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The block diagram for the on-line nonlinear algorithm implementation is shown 
in Figure 6.2. Due to the tilt coordination limit of 5 deg/sec that is needed for responses 
to surge inputs, a separate set of state equations as given in Eq. (6.26) and Riccati solver 
for Eq. (6.21) are needed for the pitch cues. 



Figure 6.2. Nonlinear Algorithm Implementation for Longitudinal Mode. 

6.2.2. Lateral Mode 

For the lateral (roll/sway) mode, the algorithm development is analogous to the 


longitudinal mode. The inputs are the same as given in Section 4.2.2 for the linear 


algorithm. The specific force / is given in Eq. (4.26). With v replacing v x , the 


differential equation given in Eq. (6.8) becomes 


v y + = 


R Sz (Aq - £>(,)«, - gu y - gA () J, u x dt + u 2 + A () u 2 


(6.28) 


which when rewritten as Eq. (6.11), will produce the state equation representation for the 
otolith model similar to Eq. (6.12), with the same system matrices except 
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D oxo = \^oro r s 0] • The state space representation for the vestibular model of the form 

of Eq. (6.13) ultimately results. Similar to Eq. (6.14), the additional state resulting from 
the optokinetic influence becomes 

Xg ~ ~T 2 Xg + T 2 ( — V 4 —Xg— GgjgTMy ) 

(6.29) 

— — T 2 x 4 — T 2 x 7 — T 2 Xg — T 2 G OTO r z U\ . 

The perceptual model representation of the same form as Eqs. (6.16) and (6.17) 
will then result. Additional motion platform states and filtered white noise states defined 
in Section 4.2.2 are again used. The remaining development is identical in form to 
Eqs.(6.18) to (6.27), resulting in state equations with a nonlinear control law and a time- 
varying feedback matrix dependent upon solution of the Riccati equation. The on-line 
implementation for this mode is identical to Figure 6.2. 

6.2.3. Vertical Mode 

For the vertical (heave) mode, the single degree-of-freedom input u = a z , with the 
specific force /j = a z - g. The otolith model given in Eq. (6.8) then becomes 

AT (6-30) 


which can then be defined in state space notation as x OTO =A OTO x OTO +B OTO u, where 
x OTO are the otolith states for this mode with x 2 = v, , and 


OTO 


-B 0 0 
1 0 


B 


OTO 


Goto (^o -®o) 


OTO 


Similar to Eq. (6.14), the additional state resulting from the optokinetic influence 
becomes x 3 =-T 2 x 3 + T 2 u . The state space perceptual model now becomes: 
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(6.31) 


with 


X OTO — ^OTO X OTO ^OTO U 
V zPE ~ ^()TO X OTO ’ 



'-*0 

0 

0 " 


Goto ( A() B () ) 

A - 

^OTO — 

1 

0 

0 

’ ^OTO — 

Goto 


0 

-T 

1 2 

-T 

1 2J 


0 


[0 1 1], 


Since this mode consists of a single translational degree-of-freedom, A v = A OT() , 
By = B oto , and C v = C OT o- This results in a perceptual model of the same form as Eq. 
(6.17). Additional motion platform states and filtered white noise states defined in the 
linear algorithm development in Eqs. (4.31) and (4.32) are again used. The remaining 
development is identical in form to Eqs. (6.18) to (6.27), resulting in state equations with 
a nonlinear control law and a time-varying feedback matrix dependent upon solution of 
the Riccati equation. The block diagram for the on-line implementation for this mode is 
shown in Figure 6.3. 



Figure 6.3. Nonlinear Algorithm Implementation for Vertical Mode. 


6.2.4. Yaw Mode 

For the yaw mode, the single degree-of-freedom input is u = y/ . The state 

equations given in Eq. (6.3) apply, with output f/> replacing 9 . The additional state 
resulting from the optokinetic influence follows a development similar to Eqs. (6.4) to 
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(6.7), with y/ PE replacing 0 PE as output. Since this mode consists of a single rotational 
input, A v = A scc , By = B scc , C v = C S(X , and D v = D scc . This results in a perceptual model 
of the same form as Eq. (6.17). Additional motion platform states and filtered white 
noise states defined in the linear algorithm development in Eqs. (4.35) and (4.36) are 
again used. The remaining development is identical in form to Eqs. (6.18) to (6.27), 
resulting in state equations with a nonlinear control law and a time-varying feedback 
matrix dependent upon solution of the Riccati equation. The block diagram for the on- 
line implementation for this mode is shown in Figure 6.4. 



Figure 6.4. Nonlinear Algorithm Implementation for Yaw Mode. 

6.3. Real Time Solution of the Riccati Equation 

Solving the nonlinear Riccati equation in Eq. (6.21) is a computational challenge 
in real time as a new solution is required at each time step. Since the solution to the 
preceding time step is available, it is advantageous to use this as an initial solution when 
computing the solution for the current time step. This would result in a more refined 
computational solution that can be produced within the real time requirement, i.e., within 
a time step. The initial Riccati equation solution to the linear optimal algorithm that is 
computed off-line in MATLAB is available and can be used as the initial solution for the 
first time step. To this end we desire a technique that assumes the initial solution is 
“close” to the computational solution at a given time step. Two techniques were 
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investigated: a Newton-Raphson method and a neurocomputing approach using a 


structured neural network. 

Blackburn [59] developed a method of solution using a Newton-Raphson 
iteration. The Riccati equation given in Eq. (6.21) is first generalized as 

g(p) = psp-pa^-a; t p-r;, ( 6 . 32 ) 

where A^ = A' + al and S = B 1 R 2 'B . If we map the square matrices G and P into 
column vectors 

g(P) = vec(G) = [G u G„ ••• Gj 

r T ( 6 . 33 ) 

p = vec(P) = [P n P 21 ••• P nn \ 


then it is shown that given an initial Riccati solution p(k) , the Newton-Raphson method 
can then be used to obtain a refined solution 

p(£ + 1) = p(£)-{ 5 % p [p(£)]} g[p(fc)], (6-34) 


where 3G/3P is the Jacobian matrix, which is shown to be 



i0(a;-sp) t + (a;-sp) t 0i , 


(6.35) 


and 0 is the Kronecker product, where A 0 B = a y B . 

Structured neural networks, first introduced by Wang and Mendel [60], are a 
special neural architecture that is customized to fit the specific matrix algebra application. 
An efficient method is employed that can take advantage of the matrix structure 
associated with the algorithm. Ham and Kostanic [61] have demonstrated this approach 
in solving a wide variety of matrix algebra problems such as matrix inversion, LU 
decomposition, and solving the algebraic matrix Lyapunov equation. 
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Ham and Collins [62] developed an approach for solving the algebraic matrix 
Riccati equation. A structured neural network is used for obtaining the computational 
solution P(t) and is shown in Figure 6.5. 



Figure 6.5. Structured Neural Network for Solving the Riccati Equation. 

The error signal \(t) in Figure 6.5 is given as 

v(t)=[p(>)SP(t)- a; t p(,)-p(,)a; -R',]z((), (6.36) 

where z (t) is an excitatory input signal. An energy function is then formulated as 

E(P)=y||v||“ (where ||v||., is the Euclidean norm of v), which is minimized using the 

method of steepest descent, resulting in a system of first-order matrix differential 
equations 

p (0 = m[K (0 v ( 0 z (0 + v ( 0 zT (0 A '« (0 - v ( 0 p t ( 0 s } ( 6 - 37 ) 

where // > 0 is the learning rate parameter, and p(/) = l'(/)z(/) as shown in Figure 6.5. 

In discrete-time form (the time step At is absorbed into the learning rate // ), the learning 
rule for each training step k becomes 

P(k + l) = P(k) + juAP(k), (6.38) 
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with the update term AP ( k ) given as 

AP(k) = [A' a (k)\(k)z(k) + \(k)z T (k)A' a (k)~ v(k)p T (k)S], (6.39) 

Ham and Collins [62] noted that even though the update term is not symmetric, 
the learning rule would still converge to the positive definite, symmetric solution. They 
noted, however, that performing an additional computation resulting in a symmetric 
update term would improve convergence: 

P(k+l) = P(k) + ^-[AP(k) + AP x (k)]. (6.40) 

Ham and Collins [62] noted that the external excitatory vector input signals z(/) 
are a set of linearly independent bi-polar vectors given as 

z (1) =[l -1 ••• -l],z (2) =[-l 1 ••• -l],z (B) =[-l -1 ••• 1], (6.41) 

where each vector z (k) is presented once to the neural network in an iteration, i.e. for one 
iteration there is a total of n presentations of the training step given in Eq. (6.40), with the 
solution P(k) updated with each training step. 

The structured neural network offers some advantages over the Newton-Raphson 
method for solutions to higher-order systems. The Newton-Raphson method requires 
matrix inversion, which resulted in singular solutions for ill-conditioned systems. Matrix 
inversion is not required for the structured neural network. Eliminating both matrix 
inversion and computation of the Jacobian matrix as a Kronecker product in turn reduces 
the computational burden. For these reasons the structured network approach is selected 
for implementation into the nonlinear algorithm. 
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6.4. Algorithm Evaluation 

The nonlinear algorithm is developed to achieve the desired motion cues at an 
update rate of 60 Hz. Since the computer image generator, which provides the out-the- 
window visual imagery to the simulator pilot, also runs at 60 Hz, the motion cues would 
be synchronous with the visual cues. However, because of the computer operating 
system, the real time operating system and the I/O system on the Langley real time 
computing system, the minimum interval must be an integer multiple of 125 psec and a 
power of 2. Therefore, a time step of 16 msec or an update rate of 62.5 Hz was selected 
for the real time implementation and the pilot tests. 

For the vertical mode based upon the integrated perception model, the off-line 
solution of the Riccati equation initially produced one closed-loop eigenvalue of zero, 
which resulted in the linear optimal control weights being very difficult to tune, resulting 
in undesirable motion cues. This eigenvalue was a result of the inclusion of the 
optokinetic channel in the heave mode algorithm formulation given in Section 6.2.3; the 
formulation based on the vestibular model alone did not produce a zero eigenvalue. 
Kalman decomposition [41] was performed on the perceptual state space model of Eq. 
(6.31), resulting in a model with one less state, with no change to the perceptual response. 
Implementation of this reduced-order model removed the uncontrollable state and in turn 
eliminated the closed-loop eigenvalue of zero. The linear optimal control weights could 
then be tuned to produce the desired specific force cue; matrices Q and R d were increased 
to produce the desired onset ramp and magnitude while the filtered white noise break 
frequency y was increased to 2()yrrad/s to eliminate false cues. 
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In Figure 6.6, the vertical mode responses for a square pulse input from the 
nonlinear algorithm are compared to the linear algorithm based upon the integrated 
perception model. A learning rate parameter // = 2 x 10" is used in computing the 
Riccati equation solution. The nonlinear weight parameters Q 2 = diag( 1 .0,2.0) result in 
the desired washout characteristics. The first diagonal term Q 2 (l,l) acts upon the 
simulator displacement, reducing the z-axis displacement. The second term Q 2 (2,2) acts 
upon the simulator velocity, reducing the offset response that follows the onset cue to an 
imperceptible magnitude. 
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Figure 6.6. Nonlinear Algorithm Vertical Mode Responses. 

A further increase of Q,(l,l) does not reduce the displacement, but reduces the 
negative cue at the end of the pulse. An additional increase of Q 2 (2,2) decreases the 
offset, but the magnitude of the negative cue remains unchanged. Figure 6.6 also shows 
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the nonlinearity a and the energy norm E(P) for the heave response. The variations 
resulting for a will in turn affect the vertical response. The energy norm reaches peak 
values at the start and end of the pulse, and approaches zero as the response is washed out 
and the platform returns to a neutral position. 

The yaw mode responses for an angular acceleration pulse doublet are shown in 
Figure 6.7. Since all closed-loop eigenvalues were nonzero, Kalman decomposition of 
the perceptual model was not required. The weight R d that acts upon the simulator yaw 
displacement was reduced to 200 to produce a faster onset cue. The nonlinear weight 
parameter Q 2 = 120 reduces the false cue at the end of the pulse. 
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Figure 6.7. Nonlinear Algorithm Yaw Mode Responses. 


In this mode an upper limit a, m , equal to 1 is placed on a that restricts the yaw 
displacement. Increasing <% nm will result in a longer sustained cue, but with an increase 
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in the yaw displacement beyond 12 degrees. Note that after a , is reached, the energy 
norm E(P) rapidly decreases, increasing again when a starts to decrease below its upper 
limit. A learning rate parameter jU = 2 x 10 6 is used in computing the Riccati equation 
solution; increasing // by one order-of-magnitude results in a discontinuity in the angular 
velocity cue when a, tm is reached. 

For the two-degree-of-freedom longitudinal mode, the initial formulation with the 
integrated perception model resulted in a higher-order system (15 th -order) that is much 
larger than either heave (6 th -order) or yaw (5 th -order). Two closed-loop eigenvalues of 
zero resulted from the off-line solution of the Riccati equation. The first originated from 
the additional simulator state 0. The second resulted from the optokinetic channel for the 
translational degree-of-freedom. Removal of the additional platform state combined with 
Kalman decomposition of the perceptual model given in Eq. (6.7) eliminates the two 
closed-loop eigenvalues of zero, reducing the system to ll th -order. 

2 

The longitudinal mode responses for a surge ramp to step input of 1 m/s 
magnitude and 3 m/s /s slope are shown in Figure 6.8. Note that the nonlinear algorithm 
produces a specific force cue very close to the linear case, but with a reduction in the x- 
axis displacement. The percent reduction in displacement compared to the linear case 
increases as a function of the aircraft surge input magnitude. The peak angular velocity is 
slightly larger than the linear case, resulting in a slight increase of the pitch angle. The 
responses for both a and the energy norm (not shown) are similar to the heave and yaw 
modes in that the nonlinearity primarily affects the cue onset, with the energy norm 
washing out to zero over time. 
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Figure 6.8. Nonlinear Algorithm Longitudinal Mode Responses to Surge Input. 

The nonlinear weight parameters Q 2 = diag(0,0.6) result in the desired washout 
characteristics. Increasing the first diagonal term Q 2 (l,l) that acts upon the x-axis 
simulator displacement from 0 to 0.1 will result in increasing magnitudes of x-axis 
displacement and tilt with additional oscillatory responses. A further increase of the 
second diagonal term Q 2 (2,2) results in an increase of the peak angular velocity. A 
learning rate parameter // = 2 x 10 6 is used in computing the solution of the Riccati 
equation; the responses are unchanged with an increase of //by one order-of-magnitude. 

Figure 6.9 shows the algorithm lateral mode responses to an aircraft half-sine 
input of 3 m/s“ peak and 5 -second duration. As with the longitudinal mode, Kalman 
decomposition was performed on the integrated perception model to eliminate one zero 
eigenvalue, and the additional simulator state (j) was removed from the algorithm 
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formulation. As with the surge response, the nonlinear algorithm produces a specific 
force cue very close to the linear case, but with a 0.1-m reduction in the y-axis 
displacement. Similar to the longitudinal mode, the peak angular velocity is slightly 
larger than the linear case, but in this case the peak roll angle is slightly smaller. 


Nonlinear Algorithm Lateral Response 
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Figure 6.9. Nonlinear Algorithm Lateral Mode Responses to Sway Input. 

A learning rate parameter // = 2 x 10" 6 was again used; again the responses are 
unchanged with an increase of // by one order-of-magnitude. The nonlinear weight 
parameters Q 2 = diag(0,0.8) produce the desired washout characteristics. The effects of 
increasing the weights are the same as with the longitudinal mode. 

The responses for a pitch doublet input of 0.1 rad/sec 2 are shown in Figure 6.10. 
A learning rate parameter // = 2 x 10 6 is used in computing the solution of the Riccati 
equation. The nonlinear weight parameter Q 2 = 1 results in a large value for a, but does 
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not produce any noticeable change in the angular velocity response. Note that the 
response is scaled (by the pitch degree-of-freedom nonlinear gain), but closely follows 
the shape of the aircraft input. 


Nonlinear Algorithm Pitch Response 
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Figure 6.10. Nonlinear Algorithm Pitch Degree-of-Freedom Responses. 

The frequency characteristics of the linear state space filter are very close to a 
unity-gain filter. Since there is no benefit from solving the Riccati equation in real time, 
the formulations shown in Eqs. (6.1) to (6.23) are replaced by unity-gain filters for both 
the pitch and roll degree-of-freedom. Figure 6.11 shows the revised implementation for 
the longitudinal and lateral modes. 
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Figure 6.11. Revised Algorithm Implementation with Unity-Gain Pitch Filter. 

The systems of first-order differential equations given for the neurocomputing 
solver in Eq. (6.40) require a numerical integration algorithm. A series of algorithms 
(Euler, 2 nd -order Adams-Bashforth, 2 nd - and 4 th -order Runge-Kutta) were evaluated. No 
improvement was noticed with the higher-order methods as compared to the Euler 
method. However, for the system state equations in Eq. (6.26), the Euler method was 
found unstable for low sampling frequencies; the 2 nd -order Runge-Kutta method resulted 
in stable results for sample rates as low as 32 Hz. 

The responses using a second neurocomputing solver developed by Wang and Wu 
[63] are sensitive to the magnitude and stiffness of the closed-loop eigenvalues, with the 
responses dependent upon the choice and structure of the activation functions. The 
approach proposed by Ham and Collins [62] utilizes a structured network without 
activation functions; the responses are more robust with respect to the closed-loop 
eigenvalues. This solver yields improved responses and convergence with less 
computational burden; only one solver iteration is required per time step. 
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During implementation on the NASA Langley real time computing system, it was 
discovered that for a time step of 16 msec, the real time requirement for the nonlinear 
algorithm was not being met. The baseline software for the nonlinear algorithm, with all 
presentations of the excitatory vector z(f) given in Eq. (6.41) for each mode, i.e. n 

presentations for an /7 th -order system, resulted in an average CPU time of 34 msec. This 
CPU time also includes the contributions of the aircraft model and control loader. 
Optimizing matrix operations by taking advantage of symmetry and sparse matrices in 
the state equations and solution of the Riccati equation resulted in an average CPU time 
of 24 msec. Reducing the number of presentations of the vector z(t) for each mode to 

the first three vectors produced a CPU time of 12 msec, which met the real time 
requirement. 

This change is reflected in the nonlinear algorithm results shown in this Section. 
The results for the yaw mode were about the same as the case of 5 presentations of z (7) . 

For the heave mode, the learning rate ju was reduced from 2 x 10 6 to 2 x 10 7 to remove a 
false cue at the end of the square pulse; the results were then very close to the case of 6 
presentations of z (7) . For the longitudinal and lateral modes the results were about the 

same as the case of 11 presentations of z (7) for each mode; the only noticeable change 

was an increase in a secondary artifact in the specific force cue after the half- sine pulse as 
shown in Figure 6.9 that remains imperceptible. 

A summary of the nonlinear algorithm parameters for each of the four modes 
(longitudinal, lateral, vertical, and yaw) is given in Appendix C. 


120 



6.5. Pilot Tuning of the Algorithm 

A computer program was developed by Telban, et al. [64] for the purpose of 
driving the NASA Langley Visual Motion Simulator (VMS) that is described in Section 
2.1. This program includes both the optimal algorithm and the nonlinear algorithm. A 
general description of the program is given along with a description and flow charts of 
each cueing algorithm. Common block variable listings and a program listing are also 
provided. Procedures for tuning the nonlinear gain coefficients are also given. 

In order to determine the nonlinear gain coefficients for each degree-of-freedom 
that resulted in the most desired pilot performance, an experienced simulator pilot 
executed a series of controlled maneuvers with the optimal and nonlinear algorithms on 
the VMS. This series of maneuvers was first executed with the polynomial gain 
coefficients determined prior to testing. Coefficients for each degree-of-freedom were 
then adjusted until the simulator pilot subjectively felt the desired perception and 
performance were reached, while ensuring that the simulator motion platform limits were 
not exceeded. 

The following maneuvers were executed for both algorithms using the nonlinear 
B757 model: 

Straight Approach and Landing (with varying wind from head to tail) 

Offset Approach and Landing (with and without turbulence) 

Takeoff from Full Stop (with and without engine failure) 

Ground Maneuvers (taxiing, effect of aircraft brakes). 

No additional tuning was needed for either the straight-in or offset approach 
maneuvers. However, both algorithms showed a tendency to exceed the actuator limits 
of the motion system with the takeoff maneuver. Reducing the surge gains for the 
optimal algorithm and both the surge and pitch gains for the nonlinear algorithm resulted 
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in platform motion within the actuator limits during the takeoff maneuvers. Table 6.1 
lists the nonlinear gains by degree-of-freedom for each algorithm. From Eq. (2.8), the 
coefficients c b c 2 , and c 3 are given for each degree-of-freedom. 


Table 6.1. Nonlinear Gain Coefficients for the Cueing Algorithms. 



Optimal Algorithm 

Non 

linear Algorithm 

Degree-of- 

Freedom 

Ci 

c 2 

c 3 

c. 

c 2 

c 3 

Surge (X) 

0.6 

-0.055 

0.002 

0.5 

-0.05 

0.002 

Sway (Y) 

0.5 

-0.055 

0.002 

0.4 

-0.035 

0.001 

Heave (Z) 

0.6 

-0.082 

0.0038 

0.6 

-0.082 

0.0038 

Roll (p) 

0.3 

-0.3 

0.1 

0.3 

-0.3 

0.1 

Pitch (q) 

0.4 

-0.54 

0.26 

0.3 

-0.3 

0.1 

Yaw (r) 

1.1 

-1.46 

0.64 

1.1 

-1.46 

0.64 


6.6. Comparison of Motion Cueing Algorithms 

Algorithm responses using test runs for each degree-of-freedom are given in 
Appendix D. Comparisons are made (with the linear optimal algorithm) of both specific 
force cues (denoted by SF in the graphs) at the pilot’s head and angular velocity cues, as 
well as the linear and angular displacement of the simulator. Actuator extension lengths 
are also compared. The number of each actuator referenced in Appendix D is shown on 
the motion platform in Figure 2.4. 

The vertical mode responses for the nonlinear algorithm for a pulse input of 1 
m/s magnitude and 10- second duration are shown in Figure 6.12. The onset ramp is 
very close to that of the adaptive and optimal algorithms, with a slightly larger peak 
magnitude. The cue is sustained for a longer duration, resulting in 33 percent more z-axis 
displacement compared to the linear optimal algorithm. The negative specific force cue 
at the end of the pulse is twice the magnitude as the adaptive algorithm response. This 
results in increased sensation that indicates the end of the pulse input. 
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Z-Axis Specific Force at Pilot Head 
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Figure 6.12. Algorithm Responses to Vertical Pulse of 1 m/s 2 Magnitude, 10-Second 
Duration. 


Figure 6.13 compares responses with the vertical pulse magnitude increased to 3 
m/s . The cue from the nonlinear algorithm is sustained for a longer duration, resulting in 
slightly less (5 percent less) z-axis displacement as compared to the linear optimal 
algorithm. The nonlinear algorithm response washes out faster due to the nonlinear 
effects generated from the Riccati equation solution. The negative cue at the end of the 
pulse is slightly smaller than the optimal algorithm response, but is much larger than the 
negative cue that results from the adaptive algorithm. 
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Z-Axis Specific Force at Pilot Head 




Time (sec) 

Figure 6.13. Algorithm Responses to Vertical Pulse of 3 m/s 2 Magnitude, 10-Second 
Duration. 

Figure 6.14 compares the algorithm responses to an aircraft longitudinal input. A 
surge ramp to step input of 1 m/s 2 peak magnitude and 3 m/s 2 /s slope is applied to each 
algorithm. The specific force response for the nonlinear algorithm does not wash out as a 
function of time, resulting from the steady-state tilt angle sustaining a constant 
magnitude. A small increase in the angular velocity (tilt) rate compared to the optimal 
algorithm is also observed. Note that the adaptive algorithm has a larger steady-state 
specific force magnitude, as well as a larger angular velocity. 
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X-Axis Specific Force at Pilot Head 
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Figure 6.14. Algorithm Responses to Surge Ramp to Step of 1 m/s 2 Magnitude, 3 
m/s 2 /s Slope. 

Figure 6.15 compares the responses from the integrated perception model for 
these surge cues. The sensed specific force, or otolith responses show the nonlinear 
algorithm closely tracks the shape of the sensed response from the aircraft. The optimal 
algorithm produces about the same onset as the nonlinear algorithm, but results in 
noticeably less sensed response, especially for the first few seconds after the peak 
magnitude is reached. The perceived velocity responses show a slightly larger magnitude 
for the nonlinear algorithm, increasing to 2 percent greater magnitude after 10 seconds. 
The adaptive algorithm shows a negative, or false specific force cue sensed at the onset 
that results in a subsequent lag and a reduction in the perceived velocity for one second. 
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X-Axis Sensed Specific Force at Pilot Head 



X-Axis Perceived Velocity at Pilot Head 
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Figure 6.15. Integrated Perception Model Responses to Surge Cues of Figure 6.14. 

Figure 6.16 compares the algorithm responses to an aircraft lateral input. A sway 
half sine of 3 m/s“ peak magnitude and 5-second duration was applied to each algorithm. 
Note that the specific force cue generated by the adaptive algorithm has some significant 
distortion. A false cue is generated at onset, resulting in a noticeable lag in the motion 
cue response. A large peak magnitude is reached, but nearly one second after the aircraft 
input reached its peak. A large residual specific force cue remains for about three 
seconds after the aircraft input ends. The response generated by the linear optimal 
algorithm shows no negative cue at the onset, a well-shaped half sine response with a less 
noticeable lag, and much less residual specific force cue. The nonlinear algorithm results 
in a peak specific force cue that is 15 percent larger than the linear optimal algorithm, 
with even less lag and almost no residual specific force cue after the half sine input ends. 
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Y-Axis Specific Force at Pilot Head 




Figure 6.16. Algorithm Responses to Sway Half Sine of 3 m/s 2 Magnitude, 5-Second 
Duration. 

Figure 6.17 compares the responses for these sway cues from the integrated 
perception model. As expected, the nonlinear algorithm peaks to a larger sensed specific 
force as compared to the optimal algorithm, resulting in a larger perceived velocity. 
After five seconds, the conflict between the vestibular and visual stimuli is reduced, 
resulting in a gradual acceptance of the visual cues governed by the optokinetic influence 
in the model. The problems noted with the adaptive algorithm are evident; the false cue 
and delayed peak are noticeable along with excessive sensed and perceived responses 
observed in the last two seconds of the pulse input. In all three algorithms, the magnitude 
of the vestibular cues eliminates the latency to onset of linearvection that would occur 
with visual stimuli alone. 


127 





if) 

E, 

>p 

'o 
_o 
CD 
> 

"D 

9 

CD 
O 

a5 

Q_ 

0123456789 10 

Time (sec) 



Figure 6.17. Integrated Perception Model Responses to Sway Cues of Figure 6.16. 

2 

The yaw mode responses for an angular acceleration doublet of 0.1 rad/s 
magnitude and 5-second duration are shown in Figure 6.18. Both the optimal and 
nonlinear algorithms extend the duration of the positive angular velocity cue about one 
second longer than the adaptive algorithm, with the nonlinear algorithm duration being 
slightly longer. Note that the false angular velocity cue near the end of the aircraft input 
is reduced for the nonlinear algorithm. The yaw angle displacement command returns to 
the neutral state (zero displacement) in less than twenty seconds, while the linear optimal 
algorithm requires more time to return to the neutral state. 
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Figure 6.18. Algorithm Responses to Yaw Doublet of 0.1 rad/s 2 and 5-Second 
Duration. 


129 




Figure 6.19 shows the roll responses for an angular acceleration doublet of 0.1 
rad/s“ magnitude and 5 -second duration. The aircraft acceleration doublet is integrated to 
produce the triangular angular velocity shown in Figure 6.19, and generates a specific 
force discontinuity at the doublet transition. Note that the specific force response for the 
nonlinear algorithm is about the same magnitude in comparison to the optimal and 
adaptive algorithm responses. 


Angular Velocity p 




Figure 6.19. Algorithm Responses to Roll Doublet of 0.1 rad/s 2 and 5-Second 
Duration. 
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Figure 6.20 shows the pitch responses of the same magnitude and duration for all 
three algorithms. Similar aircraft inputs to the roll degree-of-freedom are observed. Note 
that for a single degree-of-freedom pitch input (e.g., pitch during straight and level 
flight), the adaptive algorithm will produce the largest angular velocity and specific force 
cues, with the nonlinear algorithm cues being less than the optimal algorithm. However, 
with the addition of longitudinal cues generated during a takeoff maneuver, the aircraft 
inertial acceleration will penalize or decrease the angular velocity gain as governed by 
the cost function of Eq. (2.11). 
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Figure 6.20. Algorithm Responses to Pitch Doublet of 0.1 rad/s 2 and 5-Second 
Duration. 
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6.7. Turbulence 


Reid and Robinson [65] first addressed the problem of producing acceptable 
motion cues to turbulent gust inputs. They noted that heave is the most critical cue in 
representing turbulence, but is also the most restricted cue when constraining motion 
within the platform geometry. To overcome this limitation, they developed an approach 
in which a second set of aircraft flight equations driven only by the turbulence inputs is 
employed. The output from this augmented channel is then added to the output from the 
primary flight equations, being driven by both turbulence and the pilot control inputs, 
before serving as input to the motion system. A similar approach to that developed by 


Reid and Robinson [65] has been implemented and is shown in Figure 6.21. 
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Figure 6.21. Optimal Algorithm Vertical Mode with Augmented Turbulence 
Channel. 


The input to the augmented channel is the z-axis component w G of the turbulence 
vector G. Reid and Robinson showed that w G is the dominant turbulence component 
needed in producing vertical acceleration due to turbulence. The secondary flight 
equations can then be represented by a transfer function H a (s). The secondary 
acceleration aj, is then scaled with a constant gain K G . Both the primary and secondary 
signals are then combined before input to the vertical motion cueing filter W 22 . 
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From a simulated aircraft test run, a system identification of aircraft vertical 
accelerations in response to turbulence was performed. The transfer function H c (s) was 
then created to not only represent the acceleration, but also incorporate some desired 
motion cueing characteristics, i.e., attenuated low-frequency content and increased high- 
frequency content. The following second-order transfer function was obtained for H G (s ): 


H c (s) = 0.1 


(2.4s + l)(2.4s + l) 
(0.4.v + 1)(0. l.s + 1) ' 


(6.42) 


For the optimal algorithm, a gain of K a equal to 0.8 was chosen to maximize the 
desired sensation of turbulence while keeping the actuator extensions within the motion 
limits. Figure 6.22 shows the optimal algorithm vertical responses due to turbulence both 
with and without the augmented channel. Note that the addition of the channel results in 
larger specific force peaks along with greater z-axis displacement. Figure 6.23 shows the 
power spectral density (PSD) of the heave acceleration cues. Note that the addition of the 
augmented turbulence channel greatly increases the PSD for low and mid-range 
frequencies up to 3 Hz. 

A similar implementation to that shown in Figure 6.21 was applied to the 
nonlinear algorithm. In this approach, the linear cueing filter W 22 was replaced with the 
nonlinear heave filter, with the gain K G set equal to 1.2. Figure 6.24 shows the nonlinear 
algorithm vertical responses due to turbulence. Note that the augmented channel results 
in larger specific forces and displacements than the optimal algorithm, with a similar 
increase in the power spectral densities as shown in Figure 6.25. 
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Figure 6.22. Optimal Algorithm Motion Cues Due to Turbulence. 
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Figure 6.23. Power Spectral Density of Optimal Algorithm Motion Cues. 
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Vertical Turbulence Response 
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Figure 6.24. Nonlinear Algorithm Motion Cues Due to Turbulence. 
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Figure 6.25. Power Spectral Density of Nonlinear Algorithm Motion Cues. 
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6.8. Summary of Results 

The inclusion of the integrated perception model to the linear algorithm 
formulation will, in general, sustain the motion cues longer compared to the formulation 
solely based on the vestibular system. However, this would result in excessive simulator 
displacements that could exceed the motion system hardware limits. The addition of the 
nonlinear control law with a time-varying feedback matrix based upon simulator motion 
allows the large displacements that result from high magnitude motion cues to be washed 
out more quickly compared to cues of lower magnitude. The neurocomputing approach 
provides an effective means of updating the solution of the Riccati equation at each time 
step. Reducing the number of sub-iterations of the presentation vector z(t) results in the 
computation meeting the real time requirement, without degradation of the quality of the 
resulting motion cues. 

The vertical mode responses from the nonlinear algorithm produce a washout of 
the motion cues that significantly reduces the z-axis displacement without requiring 
additional scaling of the simulated aircraft inputs. For the longitudinal mode response to 
a surge input, the nonlinear algorithm does not produce any difference in the specific 
force cue, but shows a reduction in the x-axis displacement; the percentage reduction of 
which compared to the linear algorithm will increase with increasing aircraft inputs. The 
specific force responses shown for a large half-sine sway input to the lateral mode are 
unchanged, but again show a significant reduction with the y-axis displacement. 

The effect of the nonlinear algorithm on the yaw mode differs from the 
translational modes as the duration of the angular velocity cue is increased and a false cue 
is decreased; both effects increasing the simulator yaw displacement. A maximum limit 
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was needed on the nonlinearity a to control this increase in displacement. For pitch and 
roll degree-of-freedom inputs, the nonlinear algorithm could not yield any change in the 
motion cues. A unity-gain filter replaced the respective state space motion cueing filter 
in the longitudinal and lateral modes. 

The responses to single degree-of-freedom aircraft inputs were compared with the 
NASA adaptive algorithm and the optimal algorithm. Results for the vertical mode show 
the nonlinear algorithm producing a motion cue with a time-varying washout, sustaining 
small cues for a longer duration and washing out larger cues more quickly compared to 
the optimal algorithm. The longitudinal mode response to a surge input results in a 
specific force response with no steady-state washout due to the addition of the integrated 
perception model in the algorithm formulation. The onset of the surge response 
eliminates the false cue that persists with the adaptive algorithm. The lateral mode 
response to a sway input reveals a motion cue without the false cue or distorted shape 
observed with the adaptive algorithm, and a larger magnitude compared to the optimal 
algorithm. Yaw mode responses reveal that the nonlinear algorithm improves the motion 
cues by reducing the magnitude of negative cues and increasing the cue duration. 

In order that takeoff maneuvers be successfully completed within the motion 
system hardware limits, pilot tuning resulted in reduction of the nonlinear gain of the 
surge degree-of-freedom. This resulted in less steady-state specific force cue compared to 
the adaptive algorithm. The pitch degree-of-freedom nonlinear gain was also reduced, 
resulting in less angular velocity cues compared to the optimal algorithm. These results 
are investigated further with pilot performance testing [66]. 
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7. Conclusions and Future Research 


7.1. Conclusions 

An improved linear optimal algorithm was developed that incorporates the latest 
research on human vestibular sensation models. Using these new models, a set of linear 
motion cueing filters were synthesized and tuned using optimal control techniques. 
Preliminary pilot tests revealed unsatisfactory perception of turbulence effects. A filter 
for the vertical mode based upon a revised otolith model resulted in a significant increase 
in the magnitude of the high-frequency gain, resulting in faster responding heave cues. A 
filter for the longitudinal mode was designed with the new otolith model and resulted in 
faster responding surge cues with a reduction in the tilt coordination rate. 

The revised otolith sensation model, derived from prior research, was formulated 
with a short time constant obtained from research with afferent responses that shows one 
order-of-magnitude reduction from past work with ocular torsion responses. The 
physiological experiments from the literature produced transfer functions with a 
fractional exponent in the lead operator. By applying fractional calculus, transient 
responses to impulse and step inputs have been derived. Comparison of the transient 
response of the revised model with these responses clearly shows that a less complex 
model can generate a response that is a reasonable approximation between responses 
from the regular and irregular units. 

An integrated model of human motion perception was developed. This model 
includes models of both vestibular and visual motion sensation and incorporates the 
nonlinear interaction between the vestibular and visual stimuli. The visual estimate of 
perceived self-motion is modeled as an optokinetic influence that filters the difference 
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between the cues through a first-order low-pass filter that represents the gradual build-up 
of self-velocity. A conflict signal estimator is used to control the optokinetic influence 
gain. Models for both rotational and translational motion were developed, producing 
responses that explain the characteristics of self-motion observed in the literature. 

A nonlinear motion cueing algorithm was developed that combines features of the 
adaptive and linear optimal algorithms, and incorporates the vestibular and integrated 
perception models. A nonlinear control law was implemented that requires the solution 
of the Riccati equation in real time. The neurocomputing approach implemented for this 
task yields responses that are robust with respect to the closed-loop eigenvalues, with less 
computational burden compared to a second neurocomputing solver and a Newton- 
Raphson implementation. 

Results for the vertical mode show the nonlinear algorithm producing a motion 
cue with a time-varying washout, sustaining small cues for a longer duration and washing 
out larger cues more quickly. The addition of the integrated perception model was shown 
to improve the response to a surge input, producing a specific force response with no 
steady-state washout. Improved cues are also observed for responses to a sway input. 
The false longitudinal and lateral cues observed with the NASA adaptive algorithm were 
absent. Yaw mode responses reveal that the nonlinear algorithm improves the motion 
cues by reducing the magnitude of negative cues closer to perceptual thresholds. The 
addition of the augmented turbulence cue to the heave mode for both the optimal and 
nonlinear algorithms increases the turbulence sensation significantly so that pilot control 
inputs are influenced. 
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7.2. Suggested Future Research 

Both the optimal and nonlinear algorithm formulations resulted in rotational pitch 
filters with frequency characteristics very close to a unity-gain filter. For large aircraft 
inputs, these filters produce strong pitch cues with large pitch angles that do not wash 
out. This can result in either the simulation being stopped due to the motion limits being 
exceeded, or the motion cues being further restrained by a limiting or braking algorithm. 
Modifying the pitch filter to provide washout capability would allow large pitch angles to 
gradually decrease to a neutral position and increase the likelihood that simulations with 
large pitch angles such as takeoff be successfully completed. The pitch filter in the 
optimal algorithm longitudinal mode formulation can be modified to produce washout by 
reducing the weight component R d (4,4) that constrains the simulator pitch angle from Eq. 
(4.17). 

For the nonlinear algorithm, the corresponding weight component R d (4,4) was 
removed from the cost function to eliminate a zero closed-loop eigenvalue, resulting in 
improved convergence of the Riccati equation neurocomputing solver. Producing 
washout capability with the pitch filter would require additional research with the 
neurocomputing solver to improve convergence with a closed-loop eigenvalue of zero, or 
an augmented approach that would address this problem separately. 

A braking algorithm developed by Wu [13] was implemented for both algorithms 
and is presented with the motion cueing program implementation [64]. This braking 
algorithm did an adequate job of restraining the simulator motion as the hardware limits 
were approached, but performed poorly in returning motion control to the cueing 
algorithm. An improved algorithm that is effective in both restraining large excursions 
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and resuming regular simulator motion would allow increased nonlinear gains and 
improve motion cueing performance. One approach that is suggested is the algorithm 
developed by McFarland [67] for NASA Ames. 

The performance of the nonlinear algorithm will improve when implemented on a 
new motion system with improved actuator extensions and bandwidth. Surge and pitch 
gains can be increased to improve pilot performance. Due to the algorithm producing 
faster washout with large motion cues, the necessity for a braking algorithm to address 
large excursions may be minimal. The addition of a pitch filter with washout would 
further improve the available motion cues. This improved nonlinear algorithm could then 
be evaluated with a large homogeneous group of test pilots using the same state of the art 
techniques [66] developed for this research. 
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Appendix A. Fractional Exponent Derivation 

Because of the fractional exponent in the transfer function of Eq. (3.22), an 
elementary solution to its response cannot be readily obtained. However, an approximate 
solution to the response can be derived through the application of fractional calculus [32]. 

By first substituting the regular unit parameters into Eq. (3.22) and then 
implementing partial fraction expansion, Eq. (3.22) becomes 


1792.056 674 Qgg 5° 188 _ 0.044538 

s + 62.5 + ‘ .v + 62.5 5 + 0.0145 


0. 188 

0.016752 — . 

5 + 0.0145 


(A.l) 


In Eq. (A.l), there are two groups of two transfer functions. Each group is related 
to either the otolith mechanics (“fast”) time constant t m or the adaptation (“slow”) time 
constant t a , with one of the two transfer functions including an exponent that represents a 
fractional derivative. For the first group, the solution to the term without the fractional 
exponent can be easily obtained by taking the inverse Laplace transform of the response: 



To derive a solution to the fractional exponent term, The inverse Laplace 
transform is first obtained by applying fractional calculus [32]: 



(A. 3) 
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V ' to T{v + k + l ) 


(A. 5) 


Eq. (A. 5) is an infinite series, which for v = 0 will reduce to the Taylor series expansion 
of the exponential function: 

(citf 


L 1 


( 0 
5 


v S ~ Cl J 


L' 


f 1 A 


yS - a J 


£,(0.«) = Z 


-Jit + 1) 


e al . 


(A. 6) 


where for a = -62.5, the solution is the same as Eq. (A. 2). 

= iy_W 

/u hr{v+ky 

where /'(f) satisfies the first-order ordinary differential equation 
with the solution 


(A.7) 


(A. 8) 


f(t) = E t (v,ci) = 


at r -AM V— 1 1 

e J e u du 


r(v) 


V>1, 


(A. 9) 


Note that in Eq. (A. 9) the integral does not exist when v - 1 is less than 0. To overcome 
this problem, we use the recursion formula 


E , ( v ’ a ) = w . A + aE < ( v + *’ a ) 


r(v + i) 
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■ + 


at 


v+l 


+ ■ 


a 2 e at 


r(v + i) r(v + 2) r(v + 2) 


■ + ■ 


at 


v+l 




2 at 

a e 


;E t {v + 2,a) 
f e~ au u v+i du . 


(A. 10) 


r(v + i) r(v + 2) r(v + 2) 

Note that the integral in Eq. (A. 10) now exists since v + 1 is greater than 0. Eq. 
(A. 10) can now be used to compute the responses of the two fractional exponent transfer 
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functions given in Eq. (A.l). For each of these transfer functions, three terms are 
computed. The first two terms are analytical functions with the third term including an 
integral that requires an approximate solution. 

To evaluate the integral in Eq. (A. 10), the integral can be rewritten as 

I(t) = e- a ' [e^ u ~ n u v+x du 

= e~ at \ [e- a,(z ~ l) t v+l z v+l dz (A. 11) 

= t v+2 e~ a 'l(t). 


where / (, t ) = [V" (l Z) z v+l dz. Note that as at ->-<», the integrand in / (/) approaches 0 
for 0 < z < 1 , and 1 for z = 1. Also, for z near z = 1, we can write 


; - H = 1 +(y+ 1 x ; - 1 )+ (l,+i y~ 1 ) A -=i:c,( ; -i) j 


7=0 


(A. 12) 


This suggests we can write the integrand in I (t) as 


/(*) = ' 


at(l-z) 


'Zc J (z-iy + r , -t 1 c l (z-iy 

\j = 0 7=0 


(A. 13) 


and therefore I (t) can be rewritten as 


h t ) = Yj C j f 1 (z - 1 )' dz + R n {t), 

7=0 


(A. 14) 


where 


r #(o = i e 


at(l-z) 


V 7=0 


dz. 


(A. 15) 


and R N (t ) < K/[-at) N+l as at — > -» , where K is a constant. The integral l(t ) in Eqn. 


(A. 1 1) can now be evaluated as 
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/(r)=rv- J£c, (z-l)'dz + o -A. 

7=0 


(A. 16) 


By taking the inverse Laplace transform of Eq. (A.l) and applying Eqs. (A. 10) 
and (A. 16) to the transfer functions with fractional exponents results in the impulse 


response h(t): 


hit) = 1792.056e~ 62 ' 5 ' + 674.058 E, (-0.188,-62.5) 

- 0.044538<?~°° 145 ' - 0.016752 E, (-0.188, -0.0145) 


(A. 17) 


The response to a unit step input will now be considered. Given a system with the 
initial conditions v = 0 and x = 0 when l = 0, and an arbitrary input u(t), we look for a 
solution in the form 


x(t) = f/*(x, T)u{r)d't 


(A. 18) 


where h(x, r) is Green’s function, i.e. the system response to an impulse input, and has 
already been derived for the regular and unit transfer function in Eq. (A. 17). If we 
consider the response to a unit step, i.e. u ( r) = 1 for t > 0, the response for a term without 

the fractional exponent is simply f e aT d T = — (e at - 1) , while the response for a term with 

J() n ' > 


the fractional exponent from Eq. (A. 10) is 


f E t (v, a) dr = E t (y + 1, a). 


(A. 19) 


and applying the recursion formula in Eq. (A. 19) results in 


l( f 

E,(v + l,a) = - E t (v, a) - . 

r (^ + 1 )J 


(A. 20) 
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Applying Eqs. (A. 18) and (A. 20) to the impulse response for the regular unit parameters 


given in Table 3.1 and combining terms results in the response to a unit step: 

x(t) = 25.601 - 28.673^ 62 ' 5 ' + 3.073 A 0014493 ' 

- 10.786E, (-0.188,-62.5) + 1.156E, (-0.188,-0.014493) 


+ 9.629 


r (v + l)' 


(A. 21) 
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Appendix B. Optimal Algorithm System Parameters and Filter 
Characteristics 


Table B.l. Optimal Algorithm System Parameters 


Parameter 

Pitch/Surge 

Roll/Sway 

Yaw 

Heave 

Semicircular Canals 

Threshold (deg/sec) 

2.0 

2.0 

1.6 


Ti (sec) 

5.73 

5.73 

5.73 


t 2 (sec) 

0.005 

0.005 

0.005 


T a (sec) 

80 

80 

80 


t l (sec) 

0.06 

0.06 

0.06 


G S cc (Threshold Units) 

28.6479 

28.6479 

35.8099 


Otolith 

2 

Threshold (m/sec“) 

0.17 

0.17 


0.28 

A 0 (sec -1 ) 

0.1 

0.1 


0.1 

B 0 (sec" 1 ) 

0.2 

0.2 


0.2 

B ! (sec' 1 ) 

62.5 

62.5 



Goto (Threshold Units) 

4.7059* Bi 

4.7059* B ! 


2.8571 

Filtered White Noise Break Frequency 

A„(l,l) 

1 

1 

1 

1 

A n (2,2) 

n 

K 



Pena 

lty Weights 

0(U) 

1 

1 

1 

1 

Q(2,2) 

10 

10 



Rd(l,l) 

8 

8 

0.1 

0.1 

Rd(2,2) 

4 

4 

300 

4 

Rd(3,3) 

1 

1 


1 

Rd (4,4) 

250 

250 



R(l,l) 

1 

1 

1 

1 

R(2,2) 

1 

1 
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Table B.2. Optimal Algorithm Filters with Young-Meiry Otolith Model. 


Poles 

W n (s) Zeros 

W 12 (s) Zeros 

W 22 (s) Zeros 

Longitudinal (Pitch/Surge) Filter 

-68.1595 

-68.7884 

-200 

-83.3944 

-7.0133 

-6.5825 

-11.0387-14. 9159i 

-13.4373 

-0. 6401-0. 9364i 

-0. 6880-0. 9932i 

-11. 0387+14. 9159i 

-1.9681 

-0. 6401+0. 9364i 

-0. 6880+0. 9932i 

-0. 3930-0. 3584i 

0 

-1. 1077-0. 1855i 

-1. 0203-0. 3224i 

-0. 3930+0. 3584i 

0 

-1. 1077+0. 1855i 

-1. 0203+0. 3224i 

-0.0468 

0 

-0.0774 

-0.0753 

0 

0 

Lateral (Roll/Sway) Filter 

-59.9780 

-60.4781 

-200 

-77.7778 

-7.9333 

-7.4642 

-11.7349-13. 9453i 

-13.6704 

-0. 6034-0. 9274i 

-0. 6391-0. 9757i 

-11. 7349+13. 9453i 

-1.9129 

-0. 6034+0. 9274i 

-0. 6391+0. 9757i 

-0. 3862-0. 3524i 

0 

-1. 0354-0. 2332i 

-0.9707-0.3 186i 

-0. 3862+0. 3524i 

0 

-1. 0354+0. 2332i 

-0. 9707+0. 3186i 

-0.0440 

0 

-0.0774 

-0.0753 

0 

0 

Vertical (Heave) Filter 

-5.6454 



-48.1817 

-0. 5503-0. 4841i 



-0.1827 

-0.5503+0.4841i 



0 

-0.1865 



0 

-0.1587 



0 

Yaw Filter 

-16.6776 

-11.8526 



-0.4070 

-0.1429 



-0.1944 

0 



-0.0316 

0 
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Table B.3. Optimal Algorithm Filters with New Otolith Model. 


Poles 

W„(s) Zeros 

W 12 (s) Zeros 

W 22 (s) Zeros 

Longitudinal (Pitch/Surge) Filter 

-13.4477 

-13.0429 

-1436 

-5.5720 

-0. 5214-0. 8914i 

-0. 5317-0. 9162i 

-3.2046 

-1.5649 

-0. 5214+0. 8914i 

-0. 5317+0. 9162i 

-0. 3868-0. 3566i 

0 

-0.903 l-0.2666i 

-0. 8718-0. 2961i 

-0. 3868+0. 3566i 

0 

-0.903 l+0.2666i 

-0. 8718+0. 2961i 

-0.0762 

0 

-0.1018 

-0.0995 

0 

0 

Lateral (Roll/Sway) Filter 

-13.4477 

-13.0429 

-1436 

-5.5720 

-0. 5214-0. 8914i 

-0. 5317-0. 9162i 

-3.2046 

-1.5649 

-0. 5214+0. 8914i 

-0. 5317+0. 9162i 

-0. 3868-0. 3566i 

0 

-0.903 l-0.2666i 

-0. 8718-0. 2961i 

-0. 3868+0. 3566i 

0 

-0.903 l+0.2666i 

-0. 8718+0. 2961i 

-0.0762 

0 

-0.1018 

-0.0995 

0 

0 

Vertical (Heave) Filter 

-0.5870-0.5 120i 



-0.1943 

-0.5870+0.5 120i 



0 

-0.1993 



0 

-0.1587 



0 

Yaw Filter 

-16.6662 

-11.7582 



-0.4429 

-0.1420 



-0.1506 

0 



-0.0183 

0 
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Appendix C. Nonlinear Algorithm System Parameters 


Table C.l. Nonlinear Algorithm System Parameters 


Parameter 

Pitch/Surge 

Roll/Sway 

Yaw 

Heave 






Semicircular Canals 

Threshold (deg/sec) 

2.0 

2.0 

1.6 


T\ (sec) 

5.73 

5.73 

5.73 


r A (sec) 

80 

80 

80 


Gscc (Threshold Units) 

28.6479 

28.6479 

35.8099 


r OK (sec) 

2.0 

2.0 

2.0 


Otolith 

2 

Threshold (m/sec“) 

0.17 

0.17 


0.28 

A 0 (sec" 1 ) 

0.1 

0.1 


0.1 

B 0 (sec -1 ) 

0.2 

0.2 


0.2 

B[ (sec' 1 ) 

62.5 

62.5 



Goto (Threshold Units) 

4.7059* Bj 

4.7059* B, 


2.8571 

Cik (sec) 

2.0 

2.0 


2.0 






Filtered White N 

oise Break Frequency 

A n ll 

1 

1 

1 

20;r 

A„22 

71 

71 



Linear Weights 

0(U) 

1 

1 

1 

1 

Q(2,2) 

300 

300 



Rd(l,l) 

8 

8 

0.1 

40 

Rd(2,2) 

4 

4 

300 

400 

Rd(3,3) 

1 

1 


40 

R(l,l) 

1 

1 

1 

200 

R(2,2) 

1 

1 



Nonlinear Parameters 


2.0 x 10' 6 

2.0 x 10' 6 

2.0 x 10' 6 

2.0 x 10' 7 

Q 2(U) 

0 

0 

1.0 

1.0 

02(2,2) 

0.6 

0.8 


2.0 

^max 

1.0 

1.0 

1.0 

0.2 
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Appendix D. Optimal and Nonlinear Algorithm Comparison Figures 


X-Axis SF at Pilot Head Angular Velocity q 



Actuator Extension 


Actuator Extension 



Figure D.l. Surge Ramp to Step, 1 m/s 2 Magnitude, 3 m/s 2 /s Slope. 
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X-Axis SF at Pilot Head Angular Velocity q 



Actuator Extension Actuator Extension 



Figure D.2. Surge Ramp to Step, 2 m/s 2 Magnitude, 3 m/s 2 /s Slope. 
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Y-Axis SF at Pilot Head Angular Velocity p 



Actuator Extension Actuator Extension 




Figure D.4. Sway Half Sine, 1 m/sec 2 Magnitude, 5-Second Duration. 
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Y-Axis SF at Pilot Head Angular Velocity p 






Actuator Extension Actuator Extension 





Figure D.4. Sway Half Sine, 3 m/sec 2 Magnitude, 5-Second Duration. 
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Z-Axis Specific Force at Pilot Head 




Actuator Extension Actuator Extension 





Figure D.5. Vertical Pulse, 1 m/sec 2 Magnitude, 10-Second Duration. 
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Specific Force (m/s 


Z-Axis Specific Force at Pilot Head 



Time (sec) 


Actuator Extension Actuator Extension 





Time (sec) Time (sec) 


Figure D.6. Vertical Pulse, 3 m/sec 2 Magnitude, 10-Second Duration. 
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5 


Y-Axis SF at Pilot Head 


Angular Velocity p 





Actuator Extension 


Actuator Extension 





Figure D.7. Roll Doublet, 0.1 rad/sec 2 Magnitude, 5-Second Duration. 
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X-Axis SF at Pilot Head Angular Velocity q 




Actuator Extension Actuator Extension 



Figure D.8. Pitch Doublet, 0.1 rad/sec 2 Magnitude, 5-Second Duration. 
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Angular Velocity r 



Actuator Extension Actuator Extension 





Figure D.9. Yaw Doublet, 0.1 rad/sec 2 Magnitude, 5-Second Duration. 
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